Homoscedastic Error Variances

For the visualization of homocedastic and hetercedastic variance, we simulate:

\[ y \sim N \left ( -1 + 2x, 1 \right ) \]

The funnel-shaped heterocedasticy is simulated using:

\[ y \sim N \left ( -1 + 2x, \left ( 0.1 + 0.3 \left ( x + 3 \right ) \right )^{2} \right ) \]

par(mfrow = c(2, 2))

# Homecedastic variance
X <- seq(-2, 2, length.out = 300)

homecedastic_random <- function (x){
  rnorm(1, mean = -1 + 2 * x,  sd = 1)[1]
}

homo_y <- sapply(X, homecedastic_random)
homo_e <- (-1 + 2 * X) - homo_y

plot(X, homo_y, xlab='', ylab='', main='homocedastic variance')
abline(-1, 2, col = "red", lwd=2)
plot(X, homo_e, xlab='', ylab='', main='homocedastic variance errors')
abline(0, 0, col = "red", lwd=2)

heterocedastic_random <- function (x){
  rnorm(1, mean = -1 + 2 * x, sd = 0.1 + 0.3 * (x + 3))[1]
}

hetero_y <- sapply(X, heterocedastic_random)
hetero_e <- (-1 + 2 * X) - hetero_y

plot(X, hetero_y, xlab='', ylab='', main='funnel-shaped heterocedastic variance')
abline(-1, 2, col = "red", lwd=2)
plot(X, hetero_e, xlab='', ylab='', main='funnel-shaped heterocedastic variance errors')
abline(0, 0, col = "red", lwd=2)

Example 3.1 Munich Rent Index—Heteroscedastic Variances

par(mfrow = c(1, 2))

munich_rent_url <- "https://www.uni-goettingen.de/de/document/download/64c29c1b1fccb142cfa8f29a942a9e05.raw/rent99.raw"

munich_rent_index <- read.table(
  url(munich_rent_url),
  header = 1,
  colClasses = c(
    "numeric", "numeric", "numeric",
    "numeric", "factor", "factor",
    "factor", "factor", "factor"
  )
)

simple_reg <- lm(rent ~ area, data = munich_rent_index)

plot(munich_rent_index$area, munich_rent_index$rent, xlab = 'area in sqm', ylab = 'net rent in Euro')
abline(simple_reg, col = "red", lwd = 2)

predicted <- predict(simple_reg)

plot(munich_rent_index$area, predicted - munich_rent_index$rent, xlab = 'area in sqm', ylab = 'residuals')
abline(0, 0, col = "red", lwd = 2)

Uncorrelated Errors

Simulated autocorrelated errors with positive correlation are generated using:

\[ y_{i} = -1 + 2 x_{i} + \epsilon_{i} \] \[ \text{Where } \epsilon_{i} = 0.9 \epsilon_{i - 1} + \mu_{i} \wedge \mu_{i} \sim N(0, 0.5^{2}) \]

par(mfrow = c(1, 2))

n <- 100
X <- seq(-3, 3, length.out = n)
mu <- rnorm(n, mean = 0, sd = 0.5)

epsilon <- mu[1]
for(i in seq(2, n)){
  epsilon <- c(epsilon, 0.9 * epsilon[i - 1] + mu[i])
}

y_pred <- -1 + 2 * X + epsilon

plot(X, y_pred, xlab = 'x', ylab = '', main = 'observations and regression line')
abline(-1, 2, col = "red", lwd = 2)

plot(epsilon, xlab = 'i', ylab = '', main = 'time series of the errors')
abline(0, 0, col = "red", lwd = 2)

Simulated autocorrelated errors with negative correlation are generated using:

\[ y_{i} = -1 + 2 x_{i} + \epsilon_{i} \] \[ \text{Where } \epsilon_{i} = - 0.9 \epsilon_{i - 1} + \mu_{i} \wedge \mu_{i} \sim N(0, 0.5^{2}) \]

par(mfrow = c(1, 2))

n <- 100
X <- seq(-3, 3, length.out = n)
mu <- rnorm(n, mean = 0, sd = 0.5)

epsilon <- mu[1]
for(i in seq(2, n)){
  epsilon <- c(epsilon, -0.9 * epsilon[i - 1] + mu[i])
}

y_pred <- -1 + 2 * X + epsilon

plot(X, y_pred, xlab = 'x', ylab = '', main = 'observations and regression line')
abline(-1, 2, col = "red", lwd = 2)

plot(epsilon, xlab = 'i', ylab = '', main = 'time series of the errors')
abline(0, 0, col = "red", lwd = 2)

Mispecified model is simulated using:

\[ y_{i} = \sin(x_{i}) + x_{i} + \epsilon_{i} \] \[ \text{Where } \epsilon_{i} \sim N(0, 0.3^{2}) \]

par(mfrow = c(2, 2))

n <- 100
X <- seq(-3, 3, length.out = n)
epsilon <- rnorm(n, mean = 0, sd = 0.3)

y <- sin(X) + X + epsilon

plot(X, y, xlab = '', ylab = '', main = 'observations and true function')
lines(X, sin(X) + X, col = "red", lwd = 2)

linear_model <- lm(y ~ X)

plot(X, y, xlab = '', ylab = '', main = 'observations and regression line')
abline(linear_model, col = "red", lwd = 2)

y_pred <- predict(linear_model, data = X)
res <- y_pred - y
plot(X, res, xlab = '', ylab = '', main = 'residuals')
abline(0, 0, col = "red", lwd = 2)

Additivity of Errors

We simulate data using:

\[ y_{i} = \exp(1 + x_{i1} - x_{i2} + \epsilon_{i}) \]

With: \[ \epsilon_{i} \sim N \left ( 0, {0.4}^{2} \right ) \]

We plot this model:

par(mfrow = c(2, 2))

n <- 50

x1 <- seq(0, 3, length.out = n)
x2 <- seq(0, 3, length.out = n)
x_grid <- expand.grid(x1 = x1, x2 = x2)
epsilon <- rnorm(n^2, mean = 0, sd = 0.4)
y <- exp(1 + x_grid$x1 - x_grid$x2 + epsilon)

plot(
  x_grid$x1, y,
  main = 'scatter plot: y versus x1',
  xlab = 'x1',
  ylab = 'y'
)
plot(
  x_grid$x2, y,
  main = 'scatter plot: y versus x2',
  xlab = 'x2',
  ylab = 'y'
)

plot(
  x_grid$x1, log(y),
  main = 'scatter plot: log(y) versus x1',
  xlab = 'x1',
  ylab = 'log(y)'
)
plot(
  x_grid$x2, log(y),
  main = 'scatter plot: log(y) versus x2',
  xlab = 'x2',
  ylab = 'log(y)'
)

Example 3.2 Supermarket Scanner Data

Data not available online.

Example 3.3 Munich Rent Index — Variable Transformation

Importing data using script:

source("import_data/munich_rent_index.R")

We do the regression model:

\[ \text{rentsqm}_{i} = \beta_{0} + \beta_{1} \cdot f(\text{area}_{i}) + \epsilon_{i} \]

With both \(f(\cdot)\):

\[ f(\text{area}_{i}) = \text{area}_{i} \]

For linear model and:

\[ f(\text{area}_{i}) = \frac{1}{\text{area}_{i}} \]

For future use we define those models:

Model 1 (\(M1\)):

\[ \mathbb{E}(\text{rentqsm}) = \beta_{0} + \beta_{1} \text{area} \]

Model 2 (\(M2\)):

\[ \mathbb{E}(\text{rentqsm}) = \beta_{0} + \beta_{1} \frac{1}{\text{area}} \]

For nonlinear relationship between \(\text{rentsqm}\) and \(\text{area}\). The plotted graph below are those two models with the left panels the linear regression without the transformation and the right panels with the inverse transformation.

par(mfrow = c(3, 2))

m1 <- lm(rentsqm ~ area, data = munich_rent_index)
summary(m1)

Call:
lm(formula = rentsqm ~ area, data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9622 -1.5737 -0.1102  1.5861  9.9674 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.46883    0.12426   76.20   <2e-16 ***
area        -0.03499    0.00174  -20.11   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.291 on 3080 degrees of freedom
Multiple R-squared:  0.1161,    Adjusted R-squared:  0.1158 
F-statistic: 404.5 on 1 and 3080 DF,  p-value: < 2.2e-16
m2 <- lm(rentsqm ~ I(1 / area), data = munich_rent_index)
summary(m2)

Call:
lm(formula = rentsqm ~ I(1/area), data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.3963 -1.5733 -0.0921  1.5824 10.1287 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.7321     0.1084   43.65   <2e-16 ***
I(1/area)   140.1783     5.9287   23.64   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.241 on 3080 degrees of freedom
Multiple R-squared:  0.1536,    Adjusted R-squared:  0.1533 
F-statistic:   559 on 1 and 3080 DF,  p-value: < 2.2e-16
plot_rentsqm_area <- function(title) {
  plot(
    munich_rent_index$area,
    munich_rent_index$rentsqm,
    main = title,
    xlab = 'area in sqm',
    ylab = 'rent per sqm'
  )
}

seq_area <- seq(min(munich_rent_index$area), max(munich_rent_index$area), length.out = 100)

plot_rentsqm_area('rent per sqm vs. area')
m1_beta <- m1$coefficients
lines(seq_area, m1_beta[1] + m1_beta[2] * seq_area, col = 'red', lwd = 2)

plot_rentsqm_area('f(area) = 1/area')
m2_beta <- m2$coefficients
lines(seq_area, m2_beta[1] + m2_beta[2] * 1 / seq_area, col = 'red', lwd = 2)

plot_residuals <- function(model) {
  plot(
    munich_rent_index$area,
    model$residuals,
    main = 'residuals',
    xlab = 'area in sqm',
    ylab = 'residuals'
  )
  abline(0, 0, col = 'red', lwd = 2)
}

plot_residuals(m1)
plot_residuals(m2)

plot_av_residuals <- function(model) {
  av_residuals <- tapply(
    model$residuals, munich_rent_index$area,
    mean
  )
  plot(
    unique(munich_rent_index$area),
    av_residuals,
    main = 'average residuals',
    xlab = 'area in sqm',
    ylab = 'residuals'
  )
  abline(0, 0, col = 'red', lwd = 2)
}

plot_av_residuals(m1)
plot_av_residuals(m2)

Example 3.4 Munich Rent Index — Polynomial Regression

For polynomial regressions we adjust two different model. Of second-order (Model 3 (\(M3\))):

\[ \mathbb{E}(\text{rentqsm}) = \beta_{0} + \beta_{1} \text{area} + \beta_{2} \text{area}^{2} \]

And third-order (Model 4 (\(M4\))):

\[ \mathbb{E}(\text{rentqsm}) = \beta_{0} + \beta_{1} \text{area} + \beta_{2} \text{area}^{2} + \beta_{3} \text{area}^{3} \]

par(mfrow = c(2, 2))

m3 <- lm(rentsqm ~ area + I(area^2), data = munich_rent_index)
summary(m3)

Call:
lm(formula = rentsqm ~ area + I(area^2), data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.2150 -1.5816 -0.0915  1.5869  9.9392 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.183e+01  2.684e-01  44.081   <2e-16 ***
area        -1.057e-01  7.351e-03 -14.380   <2e-16 ***
I(area^2)    4.707e-04  4.758e-05   9.892   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.255 on 3079 degrees of freedom
Multiple R-squared:  0.1433,    Adjusted R-squared:  0.1428 
F-statistic: 257.6 on 2 and 3079 DF,  p-value: < 2.2e-16
m4 <- lm(rentsqm ~ area + I(area^2) + I(area^3), data = munich_rent_index)
summary(m4)

Call:
lm(formula = rentsqm ~ area + I(area^2) + I(area^3), data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.4269 -1.5854 -0.1101  1.5948 10.0670 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.428e+01  5.730e-01  24.915  < 2e-16 ***
area        -2.175e-01  2.432e-02  -8.946  < 2e-16 ***
I(area^2)    1.981e-03  3.167e-04   6.254 4.54e-10 ***
I(area^3)   -6.105e-06  1.266e-06  -4.823 1.49e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.247 on 3078 degrees of freedom
Multiple R-squared:  0.1497,    Adjusted R-squared:  0.1489 
F-statistic: 180.7 on 3 and 3078 DF,  p-value: < 2.2e-16
plot_rentsqm_area('rent per sqm vs. area')
m3_beta <- m3$coefficients
lines(seq_area, m3_beta[1] + m3_beta[2] * seq_area + m3_beta[3] * seq_area ^ 2, col = 'red', lwd = 2)

plot_rentsqm_area('f(area) = 1/area')
m4_beta <- m4$coefficients
lines(seq_area, m4_beta[1] + m4_beta[2] * seq_area + m4_beta[3] * seq_area ^ 2 + m4_beta[4] * seq_area ^ 3, col = 'red', lwd = 2)

plot_residuals(m3)
plot_residuals(m4)

Example 3.5 Munich Rent Index — Additive Models

We define the following aditive models. Additive model 1:

\[ \mathcal{E}(\text{rentsqm}_{i}) = \beta_{0} + \beta_{1} \cdot \frac{1}{\text{area}_{i}} + \beta_{2} \cdot \text{yearc}_{i} \]

And a second one with \(\text{yearc}\) as a polinomial of degree 3:

\[ \mathcal{E}(\text{rentsqm}_{i}) = \beta_{0} + \beta_{1} \cdot \frac{1}{\text{area}_{i}} + \beta_{2} \cdot \text{yearc}_{i} + \beta_{3} \cdot {\text{yearc}_{i}}^{2} + \beta_{4} \cdot {\text{yearc}_{i}}^{3} \]

additive_m1 <- lm(rentsqm ~ I(1/area) + yearc, data = munich_rent_index)
summary(additive_m1)

Call:
lm(formula = rentsqm ~ I(1/area) + yearc, data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.2682 -1.4894 -0.1401  1.3935  9.0582 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -65.406008   3.351088  -19.52   <2e-16 ***
I(1/area)   119.360735   5.636182   21.18   <2e-16 ***
yearc         0.036033   0.001721   20.94   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.097 on 3079 degrees of freedom
Multiple R-squared:  0.2591,    Adjusted R-squared:  0.2586 
F-statistic: 538.5 on 2 and 3079 DF,  p-value: < 2.2e-16
additive_m2 <- lm(rentsqm ~ I(1/area) + yearc + I(yearc ^ 2)+ I(yearc ^ 3), data = munich_rent_index)
summary(additive_m2)

Call:
lm(formula = rentsqm ~ I(1/area) + yearc + I(yearc^2) + I(yearc^3), 
    data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9920 -1.3705 -0.1354  1.3711  8.2834 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.942e+04  2.993e+04   0.983    0.326    
I(1/area)    1.296e+02  5.536e+00  23.406   <2e-16 ***
yearc       -4.330e+01  4.592e+01  -0.943    0.346    
I(yearc^2)   2.119e-02  2.348e-02   0.902    0.367    
I(yearc^3)  -3.445e-06  4.002e-06  -0.861    0.389    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.039 on 3077 degrees of freedom
Multiple R-squared:    0.3, Adjusted R-squared:  0.2991 
F-statistic: 329.7 on 4 and 3077 DF,  p-value: < 2.2e-16

A third model is specified using the truncated year of construction (\(\text{yearc} - 1900\)):

munich_rent_index$yearco <- munich_rent_index$yearc - 1900
additive_m3 <- lm(rentsqm ~ I(1/area) + yearco + I(yearco ^ 2)+ I(yearco ^ 3), data = munich_rent_index)
summary(additive_m3)

Call:
lm(formula = rentsqm ~ I(1/area) + yearco + I(yearco^2) + I(yearco^3), 
    data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9920 -1.3705 -0.1354  1.3711  8.2834 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.417e+00  4.779e-01  11.335  < 2e-16 ***
I(1/area)    1.296e+02  5.536e+00  23.406  < 2e-16 ***
yearco      -9.434e-02  3.384e-02  -2.788  0.00534 ** 
I(yearco^2)  1.553e-03  6.728e-04   2.308  0.02105 *  
I(yearco^3) -3.445e-06  4.002e-06  -0.861  0.38947    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.039 on 3077 degrees of freedom
Multiple R-squared:    0.3, Adjusted R-squared:  0.2991 
F-statistic: 329.7 on 4 and 3077 DF,  p-value: < 2.2e-16

As in the book, we show the effect of year of construction in the polynomial of degree 3:

par(mfrow = c(1, 2))
beta_add_m2 <- additive_m2$coefficients

yearc_seq <- seq(min(munich_rent_index$yearc), max(munich_rent_index$yearc), length.out = 100)

plot_effect <- function(beta, year_seq, title){
  plot(
    year_seq,
    beta[3] * year_seq + beta[4] * year_seq^2 + beta[5] * year_seq^3,
    main = title,
    xlab = 'year of construction',
    ylab = 'effect of year of construction',
    type = 'l',
    col = 'darkblue',
    lwd = 2
  )
}

plot_effect(beta_add_m2, yearc_seq, 'effect of year of construction')

beta_add_m3 <- additive_m3$coefficients
truncate_yearc_seq <- seq(min(munich_rent_index$yearco), max(munich_rent_index$yearco), length.out = 100)
plot_effect(beta_add_m3, truncate_yearc_seq, 'effect of year of construction, coded from 18 to 96')

A new model is defined using the orthogonal polynomial of year of construction:

munich_rent_index$areainv <- (1 / munich_rent_index$area) - mean(1 / munich_rent_index$area)
poly_year <- poly(munich_rent_index$yearco - mean(munich_rent_index$yearco), 3)
munich_rent_index$yearcc <- poly_year[, 1]
munich_rent_index$yearcc2 <- poly_year[, 2]
munich_rent_index$yearcc3 <- poly_year[, 3]
additive_m4 <- lm(rentsqm ~ areainv + yearcc + yearcc2 + yearcc3, data = munich_rent_index)
summary(additive_m4)

Call:
lm(formula = rentsqm ~ areainv + yearcc + yearcc2 + yearcc3, 
    data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9920 -1.3705 -0.1354  1.3711  8.2834 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.11126    0.03674 193.575   <2e-16 ***
areainv     129.57166    5.53582  23.406   <2e-16 ***
yearcc       43.93838    2.07260  21.200   <2e-16 ***
yearcc2      27.53892    2.05958  13.371   <2e-16 ***
yearcc3      -1.75582    2.03998  -0.861    0.389    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.039 on 3077 degrees of freedom
Multiple R-squared:    0.3, Adjusted R-squared:  0.2991 
F-statistic: 329.7 on 4 and 3077 DF,  p-value: < 2.2e-16

Now we calculate the partial residuals for \(\text{area}\) and \(\text{yearc}\) define by:

\[ \hat{\epsilon}_{\text{area}, i} = \text{rentsqm}_{i} - \hat{\beta_{0}} - \hat{f}_{2}(\text{yearc}_{i}) \]

With the effect \(f(\text{yearc}_{i})\):

\[ \hat{f}_{2}(\text{yearc}) = \beta_{2} \cdot \text{yearc}_{i} + \beta_{3} \cdot {\text{yearc}_{i}}^{2} + \beta_{4} \cdot {\text{yearc}_{i}}^{3} \]

And:

\[ \hat{\epsilon}_{\text{yearc}, i} = \text{rentsqm}_{i} - \hat{\beta_{0}} - \hat{f}_{1}(\text{area}_{i}) \]

With the effect \(f(\text{yearc}_{i})\):

\[ \hat{f}_{1}(\text{yearc}) = \beta_{1} \cdot \frac{1}{\text{area}_{i}} \]

par(mfrow = c(1, 2))

beta_add_m4 <- additive_m4$coefficients

area_seq <- seq(min(munich_rent_index$area), max(munich_rent_index$area), length.out = 100)

area_effect <- function(area){
  inv_area <- 1 / area
  cent_area <- inv_area - mean(inv_area)
  beta_add_m4[2] * cent_area
}

yearc_effect <- function(yearc){
  yearc_center <- yearc - mean(yearc)
  poly_yearc <- poly(yearc_center, 3)
  beta_add_m4[3] * poly_yearc[, 1] + beta_add_m4[4] * poly_yearc[, 2] + beta_add_m4[5] * poly_yearc[, 3]
}

residual_area <- munich_rent_index$rentsqm - beta_add_m4[1] - yearc_effect(munich_rent_index$yearc)
plot(
  munich_rent_index$area,
  residual_area,
  main = 'effect of area',
  xlab = 'area in sqm',
  ylab = 'effect of area / partial residuals'
)
lines(area_seq, area_effect(area_seq), col = 'red', lwd = 2)

residual_year <- munich_rent_index$rentsqm - beta_add_m4[1] - area_effect(munich_rent_index$area)
plot(
  munich_rent_index$yearc,
  residual_year,
  main = 'effect of year of construction',
  xlab = 'year of construction',
  ylab = 'effect of year of construction / partial residuals'
)
lines(yearc_seq, yearc_effect(yearc_seq), col = 'red', lwd = 2)

Example 3.6 Munich Rent Index — Effect Coding

We adjust the regression model with the coding values of location:

\[ \mathcc{E}(\text{rentsqm}_{i}) = \beta_{0} + \beta_{1} \cdot \text{glocation}_{1} + \beta_{2} \cdot \text{tlocation}_{1} \]

munich_rent_index$glocation <- 0
munich_rent_index$glocation[munich_rent_index$location == 2] <- 1
munich_rent_index$glocation[munich_rent_index$location == 1] <- -1

munich_rent_index$tlocation <- 0
munich_rent_index$tlocation[munich_rent_index$location == 3] <- 1
munich_rent_index$tlocation[munich_rent_index$location == 1] <- -1

cod_model <- lm(
  rentsqm ~ glocation + tlocation,
  data = munich_rent_index
)
summary(cod_model)

Call:
lm(formula = rentsqm ~ glocation + tlocation, data = munich_rent_index)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.8564 -1.8376 -0.1074  1.7157 10.4494 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.46704    0.09638  77.477  < 2e-16 ***
glocation   -0.19479    0.10445  -1.865 0.062285 .  
tlocation    0.70529    0.18558   3.800 0.000147 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.426 on 3079 degrees of freedom
Multiple R-squared:  0.008867,  Adjusted R-squared:  0.008223 
F-statistic: 13.77 on 2 and 3079 DF,  p-value: 1.109e-06

Example 3.7 Munich Rent Index — Interaction with Quality of Kitchen

We import the new updated model of 2001, available at official page of the dataset as all the datasets:

munich_rent_url_01 <- "https://www.uni-goettingen.de/de/document/download/e00d039e9c1f28ab404b27b0b14931f1.raw/rent98_00.raw"

munich_rent_index_01 <- read.table(
  url(munich_rent_url_01),
  header = 1,
  colClasses = c(
    "numeric", "numeric", "numeric",
    "numeric", "integer", "integer",
    "integer", "integer", "integer",
    "integer", "numeric", "numeric",
    "numeric"
  )
)

summary(munich_rent_index_01)
   rent_euro        rentsqm_euro          area            yearc     
 Min.   :  12.83   Min.   : 0.1841   Min.   : 20.00   Min.   :1918  
 1st Qu.: 320.86   1st Qu.: 5.2601   1st Qu.: 52.00   1st Qu.:1939  
 Median : 424.12   Median : 6.8614   Median : 66.00   Median :1960  
 Mean   : 456.94   Mean   : 7.0274   Mean   : 67.92   Mean   :1956  
 3rd Qu.: 552.40   3rd Qu.: 8.7292   3rd Qu.: 82.00   3rd Qu.:1972  
 Max.   :1837.89   Max.   :17.6688   Max.   :243.00   Max.   :1999  
   glocation        tlocation          nkitchen          pkitchen      
 Min.   :0.0000   Min.   :0.00000   Min.   :0.00000   Min.   :0.00000  
 1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000  
 Median :0.0000   Median :0.00000   Median :0.00000   Median :0.00000  
 Mean   :0.3907   Mean   :0.02516   Mean   :0.09888   Mean   :0.04004  
 3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000  
 Max.   :1.0000   Max.   :1.00000   Max.   :1.00000   Max.   :1.00000  
     eboden           year01           yearc2            yearc3         
 Min.   :0.0000   Min.   :0.0000   Min.   :3678724   Min.   :7.060e+09  
 1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:3759721   1st Qu.:7.290e+09  
 Median :0.0000   Median :0.0000   Median :3841600   Median :7.530e+09  
 Mean   :0.2986   Mean   :0.3257   Mean   :3828362   Mean   :7.493e+09  
 3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:3888784   3rd Qu.:7.670e+09  
 Max.   :1.0000   Max.   :1.0000   Max.   :3996001   Max.   :7.990e+09  
    invarea        
 Min.   :0.004115  
 1st Qu.:0.012195  
 Median :0.015152  
 Mean   :0.016904  
 3rd Qu.:0.019231  
 Max.   :0.050000  

And we adjust the model:

\[ \begin{matrix} \widehat{\text{rentsqm}}_{i} = & \beta_{0} + \beta_{1} \text{areainvc}_{i} + \beta_{2} \text{yearco}_{i} + \beta_{3} \text{yearco2}_{i} + \beta_{4} \text{yearco3}_{i} \\ & + \beta_{5} \text{year01}_{i} + \beta_{6} \text{nkitchen}_{i} + \beta_{7} \text{pkitchen}_{i} \\ & + \beta_{8} \text{nkitchen}_{i} \cdot \text{year01}_{i} + \beta_{9} \text{pkitchen}_{i} \cdot \text{year01}_{i} \end{matrix} \]

munich_rent_index_01$areainvc <- munich_rent_index_01$invarea - mean(munich_rent_index_01$invarea)
munich_rent_index_01$yearco <- munich_rent_index_01$yearc - mean(munich_rent_index_01$yearc)
poly_year_01 <- poly(munich_rent_index_01$yearco, 3)
munich_rent_index_01$yearco <- poly_year_01[, 1]
munich_rent_index_01$yearco2 <- poly_year_01[, 2]
munich_rent_index_01$yearco3 <- poly_year_01[, 3]

interaction_model <- lm(
  rentsqm_euro ~ areainvc + yearco + yearco2 + yearco3 + year01 + nkitchen + pkitchen + nkitchen * year01 + pkitchen * year01,
  data = munich_rent_index_01
)

summary(interaction_model)

Call:
lm(formula = rentsqm_euro ~ areainvc + yearco + yearco2 + yearco3 + 
    year01 + nkitchen + pkitchen + nkitchen * year01 + pkitchen * 
    year01, data = munich_rent_index_01)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.6303 -1.2647 -0.1072  1.2414 10.5233 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       6.95424    0.03841 181.072  < 2e-16 ***
areainvc        123.71308    4.42630  27.950  < 2e-16 ***
yearco           49.42282    2.02589  24.396  < 2e-16 ***
yearco2          29.45781    2.02386  14.555  < 2e-16 ***
yearco3          -1.03886    1.97176  -0.527 0.598311    
year01           -0.25174    0.06678  -3.770 0.000166 ***
nkitchen          0.91567    0.12172   7.523 6.43e-14 ***
pkitchen          1.09081    0.17849   6.111 1.07e-09 ***
year01:nkitchen   0.39826    0.20922   1.904 0.057033 .  
year01:pkitchen   0.73511    0.32864   2.237 0.025346 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.963 on 4561 degrees of freedom
Multiple R-squared:  0.3397,    Adjusted R-squared:  0.3384 
F-statistic: 260.7 on 9 and 4561 DF,  p-value: < 2.2e-16

Example 3.8 Munich Rent Index — Interaction Between Living Area and Location

The model to adjust ius defined as:

\[ \mathbb{E}(\text{rentsqm}_{i}) = \beta_{0} + \beta_{1} \cdot \text{tlocation} + \beta_{2} \cdot \overline{\frac{1}{\text{area}}} + \beta{3} \cdot \text{area} \cdot \text{tlocation} \]

We use the dataset for average or top location only and change the dummy variable \(\text{tlocation}\) to \(0\) (average) and \(1\) (top). Finally, we visualize the effect as described in the book.

par(mfrow = c(1, 2))
at_rent <- subset(munich_rent_index, location == 1 | location == 3)
at_rent$tlocation[at_rent$location == 1] <- 0
summary(at_rent)
      rent            rentsqm            area            yearc      location
 Min.   :  54.92   Min.   : 1.417   Min.   : 20.00   Min.   :1918   1:1794  
 1st Qu.: 314.62   1st Qu.: 5.276   1st Qu.: 50.00   1st Qu.:1954   2:   0  
 Median : 418.00   Median : 6.849   Median : 65.00   Median :1962   3:  78  
 Mean   : 444.22   Mean   : 7.007   Mean   : 66.02   Mean   :1959           
 3rd Qu.: 537.42   3rd Qu.: 8.652   3rd Qu.: 80.00   3rd Qu.:1972           
 Max.   :1459.23   Max.   :15.428   Max.   :155.00   Max.   :1997           
                                                                            
    bath          kitchen         cheating          district        yearco     
 Mode :logical   Mode :logical   Mode :logical   623    :  53   Min.   :18.00  
 FALSE:1761      FALSE:1785      FALSE:148       1013   :  37   1st Qu.:54.00  
 TRUE :111       TRUE :87        TRUE :1724      713    :  32   Median :62.00  
                                                 1132   :  31   Mean   :59.36  
                                                 1712   :  30   3rd Qu.:72.00  
                                                 2024   :  30   Max.   :97.00  
                                                 (Other):1659                  
    areainv               yearcc             yearcc2          
 Min.   :-0.0105206   Min.   :-0.030935   Min.   :-0.0183373  
 1st Qu.:-0.0044722   1st Qu.:-0.001862   1st Qu.:-0.0147283  
 Median :-0.0015876   Median : 0.004598   Median :-0.0067317  
 Mean   : 0.0002152   Mean   : 0.002465   Mean   :-0.0008428  
 3rd Qu.: 0.0030278   3rd Qu.: 0.012674   3rd Qu.: 0.0116794  
 Max.   : 0.0330278   Max.   : 0.032863   Max.   : 0.0537875  
                                                              
    yearcc3             glocation         tlocation      
 Min.   :-0.0186291   Min.   :-1.0000   Min.   :0.00000  
 1st Qu.:-0.0172309   1st Qu.:-1.0000   1st Qu.:0.00000  
 Median :-0.0056305   Median :-1.0000   Median :0.00000  
 Mean   :-0.0004071   Mean   :-0.9583   Mean   :0.04167  
 3rd Qu.: 0.0096124   3rd Qu.:-1.0000   3rd Qu.:0.00000  
 Max.   : 0.0588009   Max.   : 0.0000   Max.   :1.00000  
                                                         
area_location_model <- lm(
  rentsqm ~ tlocation + areainv + area:tlocation,
  data = at_rent
)
summary(area_location_model)

Call:
lm(formula = rentsqm ~ tlocation + areainv + area:tlocation, 
    data = at_rent)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.2380 -1.4435 -0.1189  1.4788  7.1162 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    6.911e+00  4.967e-02 139.131   <2e-16 ***
tlocation      7.722e-01  7.280e-01   1.061    0.289    
areainv        1.431e+02  7.434e+00  19.252   <2e-16 ***
tlocation:area 1.001e-02  8.612e-03   1.163    0.245    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.101 on 1868 degrees of freedom
Multiple R-squared:  0.1775,    Adjusted R-squared:  0.1762 
F-statistic: 134.4 on 3 and 1868 DF,  p-value: < 2.2e-16
beta_al <- area_location_model$coefficients

f1_area <- function(area){
  areainvc <- (1 / area) - mean(1 / area)
  beta_al[3] * areainvc
}

f2_area <- function(area){
  beta_al[4] * area
}

small_effect <- f1_area(area_seq)
big_effect <- beta_al[2] + f1_area(area_seq) + f2_area(area_seq)

ylim <- c(min(c(small_effect, big_effect)), max(c(small_effect, big_effect)))

plot(
  area_seq,
  small_effect,
  main = 'area effects based on a linear interaction',
  xlab = 'area in sqm',
  ylab = 'area effects',
  ylim = ylim,
  type = 'l',
  col = 'darkblue',
  lwd = 2
)
lines(area_seq, big_effect, col = 'red', lwd = 2)

plot(
  area_seq,
  beta_al[2] + f2_area(area_seq),
  main = 'varying effect of top location',
  xlab = 'area in sqm',
  ylab = 'Varying effect of top location',
  type = 'l',
  col = 'darkblue',
  lwd = 2
)

Example 3.9 Orthogonal Design

Just for exemplify, we define the orthogonal process describe as an R function:

\[ \tilde{x}^{j} = x^{j} - \widetilde{X}_{j} \left ( {\widetilde{X}_{j}}^{T} \widetilde{X}_{j} \right )^{-1} {\widetilde{X}_{j}}^{T} x^{j} \]

orthogonal_design <- function(x){
  x_out <- cbind(x[, 1] - mean(x[, 1]), matrix(0, nrow = nrow(x), ncol = ncol(x) - 1))
  for (i in 2:ncol(x)){
    x_hat <- as.matrix(x_out[, 1:(i-1)])
    x_inv <- solve(t(x_hat) %*% x_hat)
    x_out[, i] <- x[, i] - x_hat %*% x_inv %*% t(x_hat) %*% x[, i]
  }
  x_out
}

We test this function on the polynomial of the \(\text{yearc}\) variable:

poly_test <- poly(munich_rent_index$yearc, 3, raw = TRUE)
orthogonal_test <- orthogonal_design(poly_test)

# And check it
c(
  orthogonal_test[, 1] %*% orthogonal_test[, 2],
  orthogonal_test[, 2] %*% orthogonal_test[, 3],
  orthogonal_test[, 1] %*% orthogonal_test[, 3]
)
[1] 7.550556e-05 8.482947e+05 2.602290e+00

Why did it not result?

Example 3.10 Munich Rent Index — Comparison of Models Using \({R}^{2}\)

We compare the models: \(M1\), \(M2\), \(M3\), and \(M4\) using its coefficient of determination \(R^{2}\)

coeff_deter <- data.frame(
  `R Squared` = c(summary(m1)$r.squared, summary(m2)$r.squared, summary(m3)$r.squared, summary(m4)$r.squared),
  row.names = c("M1", "M2", "M3", "M4")
)
coeff_deter

Example 3.11 Simple Linear Regression

For the model:

\[ y_{i} = \beta x_{i} + \epsilon_{i} \]

We can verify:

\[ \left ( {X_{n}}^{T} X_{n} \right )^{-1} \to 0 \]

And:

\[ \max_{i=1,...,n} {x_{i}}^{T} \left ( {X_{n}}^{T} X_{n} \right )^{-1} x_{i} \to 0 \text{ for } n \to \infty \]

\(x_{i} = i\):

\[ \left ( {X_{n}}^{T} X_{n} \right )^{-1} = \left ( (1, 2, 3, \cdots, n) \left ( \begin{matrix} 1 \\ 2 \\ \vdots \\ n \end{matrix} \right ) \right )^{-1} = \left ( 1 + 4 + \cdots + n^2 \right )^{-1} \to 0 \]

And:

\[ \max_{i=1,...,n} ( 1, 2, \cdots, n ) \left ( 1 + 4 + \cdots + n^2 \right )^{-1} \left ( \begin{matrix} 1 \\ 2 \\ \vdots \\ n \end{matrix} \right ) \to 0 \text{ for } n \to \infty \]

\(x_{i} = \frac{1}{i}\):

\[ \left ( {X_{n}}^{T} X_{n} \right )^{-1} = \left ( \left (1, 0.5, 0.33, \cdots, \frac{1}{n} \right ) \left ( \begin{matrix} 1 \\ 0.5 \\ \vdots \\ \frac{1}{n} \end{matrix} \right ) \right )^{-1} = \left ( \sum_{i=1}^{n} \frac{1}{i^{2}} \right )^{-1} \not{\to} 0 \]

And:

\[ \max_{i=1,...,n} \left ( 1, 0.5, \cdots, \frac{1}{n} \right ) \left ( \sum_{i = 1}^{n} \frac{1}{i^{2}} \right )^{-1} \left ( \begin{matrix} 1 \\ 0.5 \\ \vdots \\ \frac{1}{n} \end{matrix} \right ) \not{\to} 0 \text{ for } n \to \infty \]

\(x_{i} = \frac{1}{\sqrt{i}}\):

\[ \left ( {X_{n}}^{T} X_{n} \right )^{-1} = \left ( \sum_{i = 1}^{n} \frac{1}{i} \right )^{-1} \to 0 \]

And:

\[ \max_{i=1,...,n} \left ( 1, \frac{1}{\sqrt{2}}, \cdots, \frac{1}{\sqrt{n}} \right ) \left ( \sum_{i = 1}^{n} \frac{1}{i} \right )^{-1} \begin{matrix} 1 \\ \frac{1}{\sqrt{2}} \\ \vdots \\ \frac{1}{\sqrt{n}} \end{matrix} \to 0 \]

Example 3.12 Munich Rent Index — Hypothesis Testing

The regression model to adjust is:

\[ \begin{matrix} \text{rentsqm}_{i} = & \beta_{0} + \beta_{1} \text{areainvc}_{i} + \beta_{2} \text{yearco}_{i} + \beta_{3} \text{yearco2}_{i} + \beta_{4} \text{yearco3}_{i} \\ & + \beta_{5} \text{nkitchen}_{i} + \beta_{6} \text{pkitchen}_{i} + \beta_{7} \text{year01}_{i} + \epsilon_{i} \end{matrix} \]

hypothesis_model <- lm(
  rentsqm_euro ~ areainvc + yearco + yearco2 + yearco3 + nkitchen + pkitchen + year01,
  data = munich_rent_index_01
)
summary(hypothesis_model)

Call:
lm(formula = rentsqm_euro ~ areainvc + yearco + yearco2 + yearco3 + 
    nkitchen + pkitchen + year01, data = munich_rent_index_01)

Residuals:
   Min     1Q Median     3Q    Max 
-7.674 -1.270 -0.113  1.242 10.479 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6.93239    0.03757 184.539  < 2e-16 ***
areainvc    123.76987    4.42908  27.945  < 2e-16 ***
yearco       49.37274    2.02573  24.373  < 2e-16 ***
yearco2      29.49985    2.02515  14.567  < 2e-16 ***
yearco3      -0.88418    1.97229  -0.448  0.65396    
nkitchen      1.04340    0.10151  10.279  < 2e-16 ***
pkitchen      1.30205    0.15226   8.552  < 2e-16 ***
year01       -0.18524    0.06213  -2.981  0.00288 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.964 on 4563 degrees of freedom
Multiple R-squared:  0.3385,    Adjusted R-squared:  0.3375 
F-statistic: 333.6 on 7 and 4563 DF,  p-value: < 2.2e-16

So we want to test the hypothesis:

\[ \begin{matrix} H_{0} : \beta_{7} = 0 & \text{against} & H_{1} : \beta_{7} \neq 0 \end{matrix} \]

About the significance of the follow-up measure.

\[ \begin{matrix} H_{0} : \left ( \begin{matrix} \beta_{5} \\ \beta_{6} \end{matrix} \right ) = \left ( \begin{matrix} 0 \\ 0 \end{matrix} \right ) & \text{against} & H_{1} : \left ( \begin{matrix} \beta_{5} \\ \beta_{6} \end{matrix} \right ) \neq \left ( \begin{matrix} 0 \\ 0 \end{matrix} \right ) \end{matrix} \]

About the significance of the kitchen variable. And a final test about the significance of differentiate between these two type of kitchens:

\[ \begin{matrix} H_{0} : \beta_{5} - \beta_{6} = 0 & \text{against} & H_{1} : \beta_{5} - \beta_{6} \neq 0 \end{matrix} \]

Example 3.13 Munich Rent Index — Standard Output and Hypothesis Tests

For the first test (\(H_{0} : \beta_{7} = 0\)), we calculate:

\[ t_{7} = \frac{\hat{\beta}_{7}}{{\widehat{\text{Var}(\hat{\beta}_{7})}}^{1/2}} \sim t_{1 - \alpha / 2} (n - p) \]

The covariance matrix is:

cov_matrix <- vcov(hypothesis_model)
cov_matrix
             (Intercept)     areainvc       yearco      yearco2      yearco3
(Intercept)  0.001411211  0.006823735  0.004991278  0.005173670 -0.002333528
areainvc     0.006823735 19.616785356 -1.294964483  1.420604414  0.122987111
yearco       0.004991278 -1.294964483  4.103591466  0.046451741 -0.006404590
yearco2      0.005173670  1.420604414  0.046451741  4.101229176  0.027782090
yearco3     -0.002333528  0.122987111 -0.006404590  0.027782090  3.889943976
nkitchen    -0.001093652 -0.091944737 -0.025946481 -0.026748293  0.006345756
pkitchen    -0.001128930  0.022989833 -0.042249319 -0.049085309 -0.016014007
year01      -0.001270640  0.004137400 -0.002253660 -0.001730023  0.007205398
                 nkitchen      pkitchen        year01
(Intercept) -1.093652e-03 -0.0011289298 -1.270640e-03
areainvc    -9.194474e-02  0.0229898326  4.137400e-03
yearco      -2.594648e-02 -0.0422493193 -2.253660e-03
yearco2     -2.674829e-02 -0.0490853093 -1.730023e-03
yearco3      6.345756e-03 -0.0160140070  7.205398e-03
nkitchen     1.030419e-02  0.0014042193  5.682974e-05
pkitchen     1.404219e-03  0.0231824417  1.902243e-04
year01       5.682974e-05  0.0001902243  3.860040e-03

We check the \(\text{Var}(\hat{\beta}_{7})\) is the last value.

var_7 <- cov_matrix[8, 8]
t_7 <- hypothesis_model$coefficients[8] / sqrt(var_7)
t_7
   year01 
-2.981496 

And this coincides with the value given in the summary.

For the second test, we compute the \(F\) value as:

\[ F = \frac{1}{r} {\hat{\beta}}^{T} \widehat{\text{Cov} \left ( \hat{\beta} \right )}^{-1} \hat{\beta} \sim F_{1, n - p} \]

With the \(\hat{\beta}\) define in the test \(\left ( \begin{matrix} \hat{\beta_{5}} \\ \hat{\beta_{6}} \end{matrix} \right )\), then \(r=2\) and \(n - p = n - 8\):

df <- nrow(munich_rent_index_01) - 8
r <- 2
beta_hat <- as.matrix(hypothesis_model$coefficients[6:7])
f_56 <- (1/2) * t(beta_hat) %*% solve(cov_matrix[6:7, 6:7]) %*% beta_hat
f_56
         [,1]
[1,] 82.08307

We get expected \(F\):

f_crit <- qf(p = 0.95, r, df)
f_crit
[1] 2.9977

And we reject the null hypothesis.

The third and final test comparing \(\beta_{5}\) and \(\beta_{6}\) between each other. We need:

\[ \widehat{\text{Var}(\hat{\beta}_{5} - \hat{\beta}_{6})} = \widehat{\text{Var}(\hat{\beta}_{5})} + \widehat{\text{Var}(\hat{\beta}_{6})} - 2 \cdot \widehat{\text{Cov}(\hat{\beta}_{5}, \hat{\beta}_{6})} \]

Using the \(\widehat{\text{Cov}(\hat{\beta})}\) matrix:

cov_5_6 <- cov_matrix[6, 7]
var_5 <- cov_matrix[6, 6]
var_6 <- cov_matrix[7, 7]
var_5_6 <- var_5 + var_6 - 2 * cov_5_6
f_5_6 <- (hypothesis_model$coefficients[6] - hypothesis_model$coefficients[7])^2 / var_5_6
f_5_6
nkitchen 
 2.18071 

And this \(F\) critical, using \(r=1\):

f_crit <- qf(p = 0.95, 1, df)
f_crit
[1] 3.843498

So, in this case, we do not reject the null hypothesis.

Example 3.14 Munich Rent Index — Confidence Intervals

The confidence interval for \(\beta_{7}\) is obtain using:

\[ \beta_{7} \pm t_{n - p} \left (1 - \alpha / 2 \right ) \cdot \widehat{\text{Var}(\hat{\beta}_{7})}^{1/2} \]

inter <- qt(1 - 0.05 / 2, df) * sqrt(var_7)
c(hypothesis_model$coefficients[8] - inter, hypothesis_model$coefficients[8] + inter)
    year01     year01 
-0.3070414 -0.0634347 

Example 3.15 Munich Rent Index — Prediction Intervals

Using the prediction interval:

\[ {x_{0}}^{T} \cdot \hat{\beta} \pm t_{n-p}(1 - \alpha / 2) \hat{\sigma} \left ( 1 + {x_{0}}^{T} \left (X^{T} X \right )^{1} x_{0} \right )^{1/2} \]

We can obtain the \(\hat{\sigma}\) directly from the model or using:

\[ (n - p) \hat{\sigma}^2 = \epsilon^{T} \epsilon \]

sigma_model <- summary(hypothesis_model)$sigma
sigma_2 <- t(as.matrix(hypothesis_model$residuals)) %*% as.matrix(hypothesis_model$residuals) / (nrow(munich_rent_index_01) - 8)
sigma <- sqrt(sigma_2)
c(sigma_model, sigma)
[1] 1.964112 1.964112

We are also going to use the design matrix:

design_matrix <- as.matrix(hypothesis_model$model)
design_matrix[, 1] <- 1
colnames(design_matrix)[1] <- "1"
head(design_matrix)
  1      areainvc       yearco      yearco2      yearco3 nkitchen pkitchen
1 1  0.0116672025 -0.011568466 -0.010460264  0.030118199        0        0
2 1 -0.0072887975 -0.011568466 -0.010460264  0.030118199        0        0
3 1  0.0175786025  0.009594692 -0.004320479 -0.013940551        0        0
4 1  0.0087368025  0.010256041 -0.003177207 -0.014569233        0        0
5 1 -0.0065948975  0.018853574  0.016932467 -0.005009151        0        0
6 1 -0.0007751975  0.003642554 -0.012015191 -0.002877678        0        0
  year01
1      0
2      0
3      0
4      0
5      0
6      0
area_seq <- seq(min(munich_rent_index_01$area), max(munich_rent_index_01$area), length.out = 100)
areainvc <- (1 / area_seq) - mean(1 / area_seq)

betam <- hypothesis_model$coefficients
yearcc <- munich_rent_index_01[munich_rent_index_01$yearc == 1918, ][1, ]
# Three last terms were added for completeness
constant <- betam[1] + betam[3] * yearcc$yearco + betam[4] * yearcc$yearco2 + betam[5] * yearcc$yearco3 + betam[6] * 0 + betam[7] * 0 + betam[8] * 0

plot(
  munich_rent_index_01$area,
  munich_rent_index_01$rentsqm_euro,
  xlab = 'area in sqm',
  ylab = 'estimated rent per sqm'
)
beta_1 <- hypothesis_model$coefficients[2]
lines(area_seq, constant + beta_1 * areainvc, col = 'red', lwd = 2)

conf_inter <- qt(1 - 0.05 / 2, df) * sqrt(cov_matrix[2, 2])
beta_conf <- c(beta_1 - conf_inter, beta_1 + conf_inter)
lines(area_seq, constant + beta_conf[1] * areainvc, col = 'darkblue', lwd = 2)
lines(area_seq, constant + beta_conf[2] * areainvc, col = 'darkblue', lwd = 2)

x_o <- as.matrix(cbind(
  1,
  areainvc,
  yearcc$yearco,
  yearcc$yearco2,
  yearcc$yearco3,
  0,
  0,
  0
))

pred_inter <- NULL

for (i in seq_len(nrow(x_o)))
  x_o_i <- as.matrix(x_o[i, ])
  pred_inter <- c(pred_inter, qt(1 - 0.05 / 2, df) * sigma * sqrt(
    1 + (
      t(x_o_i) %*% solve(t(design_matrix)%*%design_matrix) %*% x_o_i
   ))
)

lines(area_seq, constant + beta_conf[1] * areainvc - pred_inter, col = 'gray', lwd = 2)
lines(area_seq, constant + beta_conf[2] * areainvc + pred_inter, col = 'gray', lwd = 2)

This can also be done using predict:

x_o <- as.data.frame(cbind(
  areainvc,
  yearcc$yearco,
  yearcc$yearco2,
  yearcc$yearco3,
  0,
  0,
  0
))
colnames(x_o) <- colnames(hypothesis_model$model)[2:8]
confidence_interval <- predict(hypothesis_model, x_o, interval = "confidence")

plot(
    munich_rent_index_01$area,
    munich_rent_index_01$rentsqm_euro,
    xlab = 'area in sqm',
    ylab = 'estimated rent per sqm'
)
lines(area_seq, confidence_interval[, "fit"], col = 'red', lwd = 2)
lines(area_seq, confidence_interval[, "lwr"], col = 'darkblue', lwd = 2)
lines(area_seq, confidence_interval[, "upr"], col = 'darkblue', lwd = 2)

prediction_interval <- predict(hypothesis_model, x_o, interval = "prediction")
lines(area_seq, prediction_interval[, "lwr"], col = 'gray', lwd = 2)
lines(area_seq, prediction_interval[, "upr"], col = 'gray', lwd = 2)

Example 3.16 Polynomial Regression

LS0tCnRpdGxlOiAiVGhlIENsYXNzaWNhbCBMaW5lYXJNb2RlbCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyMgSG9tb3NjZWRhc3RpYyBFcnJvciBWYXJpYW5jZXMKCkZvciB0aGUgdmlzdWFsaXphdGlvbiBvZiBob21vY2VkYXN0aWMgYW5kIGhldGVyY2VkYXN0aWMgdmFyaWFuY2UsIHdlIHNpbXVsYXRlOgoKJCQKeSBcc2ltIE4gXGxlZnQgKCAtMSArIDJ4LCAxIFxyaWdodCApCiQkCgpUaGUgZnVubmVsLXNoYXBlZCBoZXRlcm9jZWRhc3RpY3kgaXMgc2ltdWxhdGVkIHVzaW5nOgoKJCQKeSBcc2ltIE4gXGxlZnQgKCAtMSArIDJ4LCBcbGVmdCAoIDAuMSArIDAuMyBcbGVmdCAoIHggKyAzIFxyaWdodCApIFxyaWdodCApXnsyfSBccmlnaHQgKQokJAoKYGBge3IsIGZpZy53aWR0aD0xMCwgZmlnLmhlaWdodD04LCBmaWcuYWxpZ249J2NlbnRlcid9CnBhcihtZnJvdyA9IGMoMiwgMikpCgojIEhvbWVjZWRhc3RpYyB2YXJpYW5jZQpYIDwtIHNlcSgtMiwgMiwgbGVuZ3RoLm91dCA9IDMwMCkKCmhvbWVjZWRhc3RpY19yYW5kb20gPC0gZnVuY3Rpb24gKHgpewogIHJub3JtKDEsIG1lYW4gPSAtMSArIDIgKiB4LCAgc2QgPSAxKVsxXQp9Cgpob21vX3kgPC0gc2FwcGx5KFgsIGhvbWVjZWRhc3RpY19yYW5kb20pCmhvbW9fZSA8LSAoLTEgKyAyICogWCkgLSBob21vX3kKCnBsb3QoWCwgaG9tb195LCB4bGFiPScnLCB5bGFiPScnLCBtYWluPSdob21vY2VkYXN0aWMgdmFyaWFuY2UnKQphYmxpbmUoLTEsIDIsIGNvbCA9ICJyZWQiLCBsd2Q9MikKcGxvdChYLCBob21vX2UsIHhsYWI9JycsIHlsYWI9JycsIG1haW49J2hvbW9jZWRhc3RpYyB2YXJpYW5jZSBlcnJvcnMnKQphYmxpbmUoMCwgMCwgY29sID0gInJlZCIsIGx3ZD0yKQoKaGV0ZXJvY2VkYXN0aWNfcmFuZG9tIDwtIGZ1bmN0aW9uICh4KXsKICBybm9ybSgxLCBtZWFuID0gLTEgKyAyICogeCwgc2QgPSAwLjEgKyAwLjMgKiAoeCArIDMpKVsxXQp9CgpoZXRlcm9feSA8LSBzYXBwbHkoWCwgaGV0ZXJvY2VkYXN0aWNfcmFuZG9tKQpoZXRlcm9fZSA8LSAoLTEgKyAyICogWCkgLSBoZXRlcm9feQoKcGxvdChYLCBoZXRlcm9feSwgeGxhYj0nJywgeWxhYj0nJywgbWFpbj0nZnVubmVsLXNoYXBlZCBoZXRlcm9jZWRhc3RpYyB2YXJpYW5jZScpCmFibGluZSgtMSwgMiwgY29sID0gInJlZCIsIGx3ZD0yKQpwbG90KFgsIGhldGVyb19lLCB4bGFiPScnLCB5bGFiPScnLCBtYWluPSdmdW5uZWwtc2hhcGVkIGhldGVyb2NlZGFzdGljIHZhcmlhbmNlIGVycm9ycycpCmFibGluZSgwLCAwLCBjb2wgPSAicmVkIiwgbHdkPTIpCmBgYAoKIyMgRXhhbXBsZSAzLjEgTXVuaWNoIFJlbnQgSW5kZXjigJRIZXRlcm9zY2VkYXN0aWMgVmFyaWFuY2VzCgpgYGB7ciwgZmlnLndpZHRoPTEwLCBmaWcuaGVpZ2h0PTQsIGZpZy5hbGlnbj0nY2VudGVyJ30KcGFyKG1mcm93ID0gYygxLCAyKSkKCm11bmljaF9yZW50X3VybCA8LSAiaHR0cHM6Ly93d3cudW5pLWdvZXR0aW5nZW4uZGUvZGUvZG9jdW1lbnQvZG93bmxvYWQvNjRjMjljMWIxZmNjYjE0MmNmYThmMjlhOTQyYTllMDUucmF3L3JlbnQ5OS5yYXciCgptdW5pY2hfcmVudF9pbmRleCA8LSByZWFkLnRhYmxlKAogIHVybChtdW5pY2hfcmVudF91cmwpLAogIGhlYWRlciA9IDEsCiAgY29sQ2xhc3NlcyA9IGMoCiAgICAibnVtZXJpYyIsICJudW1lcmljIiwgIm51bWVyaWMiLAogICAgIm51bWVyaWMiLCAiZmFjdG9yIiwgImZhY3RvciIsCiAgICAiZmFjdG9yIiwgImZhY3RvciIsICJmYWN0b3IiCiAgKQopCgpzaW1wbGVfcmVnIDwtIGxtKHJlbnQgfiBhcmVhLCBkYXRhID0gbXVuaWNoX3JlbnRfaW5kZXgpCgpwbG90KG11bmljaF9yZW50X2luZGV4JGFyZWEsIG11bmljaF9yZW50X2luZGV4JHJlbnQsIHhsYWIgPSAnYXJlYSBpbiBzcW0nLCB5bGFiID0gJ25ldCByZW50IGluIEV1cm8nKQphYmxpbmUoc2ltcGxlX3JlZywgY29sID0gInJlZCIsIGx3ZCA9IDIpCgpwcmVkaWN0ZWQgPC0gcHJlZGljdChzaW1wbGVfcmVnKQoKcGxvdChtdW5pY2hfcmVudF9pbmRleCRhcmVhLCBwcmVkaWN0ZWQgLSBtdW5pY2hfcmVudF9pbmRleCRyZW50LCB4bGFiID0gJ2FyZWEgaW4gc3FtJywgeWxhYiA9ICdyZXNpZHVhbHMnKQphYmxpbmUoMCwgMCwgY29sID0gInJlZCIsIGx3ZCA9IDIpCmBgYAoKIyMgVW5jb3JyZWxhdGVkIEVycm9ycwoKU2ltdWxhdGVkIGF1dG9jb3JyZWxhdGVkIGVycm9ycyB3aXRoIHBvc2l0aXZlIGNvcnJlbGF0aW9uIGFyZSBnZW5lcmF0ZWQgdXNpbmc6CgokJAp5X3tpfSA9IC0xICsgMiB4X3tpfSArIFxlcHNpbG9uX3tpfQokJAokJApcdGV4dHtXaGVyZSB9IFxlcHNpbG9uX3tpfSA9IDAuOSBcZXBzaWxvbl97aSAtIDF9ICsgXG11X3tpfSBcd2VkZ2UgXG11X3tpfSBcc2ltIE4oMCwgMC41XnsyfSkKJCQKCmBgYHtyLCBmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9NCwgZmlnLmFsaWduPSdjZW50ZXInfQpwYXIobWZyb3cgPSBjKDEsIDIpKQoKbiA8LSAxMDAKWCA8LSBzZXEoLTMsIDMsIGxlbmd0aC5vdXQgPSBuKQptdSA8LSBybm9ybShuLCBtZWFuID0gMCwgc2QgPSAwLjUpCgplcHNpbG9uIDwtIG11WzFdCmZvcihpIGluIHNlcSgyLCBuKSl7CiAgZXBzaWxvbiA8LSBjKGVwc2lsb24sIDAuOSAqIGVwc2lsb25baSAtIDFdICsgbXVbaV0pCn0KCnlfcHJlZCA8LSAtMSArIDIgKiBYICsgZXBzaWxvbgoKcGxvdChYLCB5X3ByZWQsIHhsYWIgPSAneCcsIHlsYWIgPSAnJywgbWFpbiA9ICdvYnNlcnZhdGlvbnMgYW5kIHJlZ3Jlc3Npb24gbGluZScpCmFibGluZSgtMSwgMiwgY29sID0gInJlZCIsIGx3ZCA9IDIpCgpwbG90KGVwc2lsb24sIHhsYWIgPSAnaScsIHlsYWIgPSAnJywgbWFpbiA9ICd0aW1lIHNlcmllcyBvZiB0aGUgZXJyb3JzJykKYWJsaW5lKDAsIDAsIGNvbCA9ICJyZWQiLCBsd2QgPSAyKQpgYGAKClNpbXVsYXRlZCBhdXRvY29ycmVsYXRlZCBlcnJvcnMgd2l0aCBuZWdhdGl2ZSBjb3JyZWxhdGlvbiBhcmUgZ2VuZXJhdGVkIHVzaW5nOgoKJCQKeV97aX0gPSAtMSArIDIgeF97aX0gKyBcZXBzaWxvbl97aX0KJCQKJCQKXHRleHR7V2hlcmUgfSBcZXBzaWxvbl97aX0gPSAtIDAuOSBcZXBzaWxvbl97aSAtIDF9ICsgXG11X3tpfSBcd2VkZ2UgXG11X3tpfSBcc2ltIE4oMCwgMC41XnsyfSkKJCQKCmBgYHtyLCBmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9NCwgZmlnLmFsaWduPSdjZW50ZXInfQpwYXIobWZyb3cgPSBjKDEsIDIpKQoKbiA8LSAxMDAKWCA8LSBzZXEoLTMsIDMsIGxlbmd0aC5vdXQgPSBuKQptdSA8LSBybm9ybShuLCBtZWFuID0gMCwgc2QgPSAwLjUpCgplcHNpbG9uIDwtIG11WzFdCmZvcihpIGluIHNlcSgyLCBuKSl7CiAgZXBzaWxvbiA8LSBjKGVwc2lsb24sIC0wLjkgKiBlcHNpbG9uW2kgLSAxXSArIG11W2ldKQp9Cgp5X3ByZWQgPC0gLTEgKyAyICogWCArIGVwc2lsb24KCnBsb3QoWCwgeV9wcmVkLCB4bGFiID0gJ3gnLCB5bGFiID0gJycsIG1haW4gPSAnb2JzZXJ2YXRpb25zIGFuZCByZWdyZXNzaW9uIGxpbmUnKQphYmxpbmUoLTEsIDIsIGNvbCA9ICJyZWQiLCBsd2QgPSAyKQoKcGxvdChlcHNpbG9uLCB4bGFiID0gJ2knLCB5bGFiID0gJycsIG1haW4gPSAndGltZSBzZXJpZXMgb2YgdGhlIGVycm9ycycpCmFibGluZSgwLCAwLCBjb2wgPSAicmVkIiwgbHdkID0gMikKYGBgCgpNaXNwZWNpZmllZCBtb2RlbCBpcyBzaW11bGF0ZWQgdXNpbmc6CgokJAp5X3tpfSA9IFxzaW4oeF97aX0pICsgeF97aX0gKyBcZXBzaWxvbl97aX0KJCQKJCQKXHRleHR7V2hlcmUgfSBcZXBzaWxvbl97aX0gXHNpbSBOKDAsIDAuM157Mn0pCiQkCgpgYGB7ciwgZmlnLndpZHRoPTEwLCBmaWcuaGVpZ2h0PTgsIGZpZy5hbGlnbj0nY2VudGVyJ30KcGFyKG1mcm93ID0gYygyLCAyKSkKCm4gPC0gMTAwClggPC0gc2VxKC0zLCAzLCBsZW5ndGgub3V0ID0gbikKZXBzaWxvbiA8LSBybm9ybShuLCBtZWFuID0gMCwgc2QgPSAwLjMpCgp5IDwtIHNpbihYKSArIFggKyBlcHNpbG9uCgpwbG90KFgsIHksIHhsYWIgPSAnJywgeWxhYiA9ICcnLCBtYWluID0gJ29ic2VydmF0aW9ucyBhbmQgdHJ1ZSBmdW5jdGlvbicpCmxpbmVzKFgsIHNpbihYKSArIFgsIGNvbCA9ICJyZWQiLCBsd2QgPSAyKQoKbGluZWFyX21vZGVsIDwtIGxtKHkgfiBYKQoKcGxvdChYLCB5LCB4bGFiID0gJycsIHlsYWIgPSAnJywgbWFpbiA9ICdvYnNlcnZhdGlvbnMgYW5kIHJlZ3Jlc3Npb24gbGluZScpCmFibGluZShsaW5lYXJfbW9kZWwsIGNvbCA9ICJyZWQiLCBsd2QgPSAyKQoKeV9wcmVkIDwtIHByZWRpY3QobGluZWFyX21vZGVsLCBkYXRhID0gWCkKcmVzIDwtIHlfcHJlZCAtIHkKcGxvdChYLCByZXMsIHhsYWIgPSAnJywgeWxhYiA9ICcnLCBtYWluID0gJ3Jlc2lkdWFscycpCmFibGluZSgwLCAwLCBjb2wgPSAicmVkIiwgbHdkID0gMikKYGBgCgojIyBBZGRpdGl2aXR5IG9mIEVycm9ycwoKV2Ugc2ltdWxhdGUgZGF0YSB1c2luZzoKCiQkCnlfe2l9ID0gXGV4cCgxICsgeF97aTF9IC0geF97aTJ9ICsgXGVwc2lsb25fe2l9KQokJAoKV2l0aDoKJCQKXGVwc2lsb25fe2l9IFxzaW0gTiBcbGVmdCAoIDAsIHswLjR9XnsyfSBccmlnaHQgKQokJAoKV2UgcGxvdCB0aGlzIG1vZGVsOgoKYGBge3IsIGZpZy53aWR0aD0xMCwgZmlnLmhlaWdodD04LCBmaWcuYWxpZ249J2NlbnRlcid9CnBhcihtZnJvdyA9IGMoMiwgMikpCgpuIDwtIDUwCgp4MSA8LSBzZXEoMCwgMywgbGVuZ3RoLm91dCA9IG4pCngyIDwtIHNlcSgwLCAzLCBsZW5ndGgub3V0ID0gbikKeF9ncmlkIDwtIGV4cGFuZC5ncmlkKHgxID0geDEsIHgyID0geDIpCmVwc2lsb24gPC0gcm5vcm0obl4yLCBtZWFuID0gMCwgc2QgPSAwLjQpCnkgPC0gZXhwKDEgKyB4X2dyaWQkeDEgLSB4X2dyaWQkeDIgKyBlcHNpbG9uKQoKcGxvdCgKICB4X2dyaWQkeDEsIHksCiAgbWFpbiA9ICdzY2F0dGVyIHBsb3Q6IHkgdmVyc3VzIHgxJywKICB4bGFiID0gJ3gxJywKICB5bGFiID0gJ3knCikKcGxvdCgKICB4X2dyaWQkeDIsIHksCiAgbWFpbiA9ICdzY2F0dGVyIHBsb3Q6IHkgdmVyc3VzIHgyJywKICB4bGFiID0gJ3gyJywKICB5bGFiID0gJ3knCikKCnBsb3QoCiAgeF9ncmlkJHgxLCBsb2coeSksCiAgbWFpbiA9ICdzY2F0dGVyIHBsb3Q6IGxvZyh5KSB2ZXJzdXMgeDEnLAogIHhsYWIgPSAneDEnLAogIHlsYWIgPSAnbG9nKHkpJwopCnBsb3QoCiAgeF9ncmlkJHgyLCBsb2coeSksCiAgbWFpbiA9ICdzY2F0dGVyIHBsb3Q6IGxvZyh5KSB2ZXJzdXMgeDInLAogIHhsYWIgPSAneDInLAogIHlsYWIgPSAnbG9nKHkpJwopCmBgYAoKIyMgRXhhbXBsZSAzLjIgU3VwZXJtYXJrZXQgU2Nhbm5lciBEYXRhCgpEYXRhIG5vdCBhdmFpbGFibGUgb25saW5lLgoKIyMgRXhhbXBsZSAzLjMgTXVuaWNoIFJlbnQgSW5kZXgg4oCUIFZhcmlhYmxlIFRyYW5zZm9ybWF0aW9uCgpJbXBvcnRpbmcgZGF0YSB1c2luZyBzY3JpcHQ6CgpgYGB7cn0Kc291cmNlKCJpbXBvcnRfZGF0YS9tdW5pY2hfcmVudF9pbmRleC5SIikKYGBgCgpXZSBkbyB0aGUgcmVncmVzc2lvbiBtb2RlbDoKCiQkClx0ZXh0e3JlbnRzcW19X3tpfSA9IFxiZXRhX3swfSArIFxiZXRhX3sxfSBcY2RvdCBmKFx0ZXh0e2FyZWF9X3tpfSkgKyBcZXBzaWxvbl97aX0KJCQKCldpdGggYm90aCAkZihcY2RvdCkkOgoKJCQKZihcdGV4dHthcmVhfV97aX0pID0gXHRleHR7YXJlYX1fe2l9CiQkCgpGb3IgbGluZWFyIG1vZGVsIGFuZDoKCiQkCmYoXHRleHR7YXJlYX1fe2l9KSA9IFxmcmFjezF9e1x0ZXh0e2FyZWF9X3tpfX0KJCQKCkZvciBmdXR1cmUgdXNlIHdlIGRlZmluZSB0aG9zZSBtb2RlbHM6CgoqKk1vZGVsIDEgKCRNMSQpKio6CgokJApcbWF0aGJie0V9KFx0ZXh0e3JlbnRxc219KSA9IFxiZXRhX3swfSArIFxiZXRhX3sxfSBcdGV4dHthcmVhfQokJAoKKipNb2RlbCAyICgkTTIkKSoqOgoKJCQKXG1hdGhiYntFfShcdGV4dHtyZW50cXNtfSkgPSBcYmV0YV97MH0gKyBcYmV0YV97MX0gXGZyYWN7MX17XHRleHR7YXJlYX19CiQkCgpGb3Igbm9ubGluZWFyIHJlbGF0aW9uc2hpcCBiZXR3ZWVuICRcdGV4dHtyZW50c3FtfSQgYW5kICRcdGV4dHthcmVhfSQuIFRoZSBwbG90dGVkIGdyYXBoIGJlbG93IGFyZSB0aG9zZSB0d28gbW9kZWxzIHdpdGggdGhlIGxlZnQgcGFuZWxzIHRoZSBsaW5lYXIgcmVncmVzc2lvbiB3aXRob3V0IHRoZSB0cmFuc2Zvcm1hdGlvbiBhbmQgdGhlIHJpZ2h0IHBhbmVscyB3aXRoIHRoZSBpbnZlcnNlIHRyYW5zZm9ybWF0aW9uLgoKYGBge3IsIGZpZy53aWR0aD0xMCwgZmlnLmhlaWdodD0xMiwgZmlnLmFsaWduPSdjZW50ZXInfQpwYXIobWZyb3cgPSBjKDMsIDIpKQoKbTEgPC0gbG0ocmVudHNxbSB+IGFyZWEsIGRhdGEgPSBtdW5pY2hfcmVudF9pbmRleCkKc3VtbWFyeShtMSkKCm0yIDwtIGxtKHJlbnRzcW0gfiBJKDEgLyBhcmVhKSwgZGF0YSA9IG11bmljaF9yZW50X2luZGV4KQpzdW1tYXJ5KG0yKQoKcGxvdF9yZW50c3FtX2FyZWEgPC0gZnVuY3Rpb24odGl0bGUpIHsKICBwbG90KAogICAgbXVuaWNoX3JlbnRfaW5kZXgkYXJlYSwKICAgIG11bmljaF9yZW50X2luZGV4JHJlbnRzcW0sCiAgICBtYWluID0gdGl0bGUsCiAgICB4bGFiID0gJ2FyZWEgaW4gc3FtJywKICAgIHlsYWIgPSAncmVudCBwZXIgc3FtJwogICkKfQoKc2VxX2FyZWEgPC0gc2VxKG1pbihtdW5pY2hfcmVudF9pbmRleCRhcmVhKSwgbWF4KG11bmljaF9yZW50X2luZGV4JGFyZWEpLCBsZW5ndGgub3V0ID0gMTAwKQoKcGxvdF9yZW50c3FtX2FyZWEoJ3JlbnQgcGVyIHNxbSB2cy4gYXJlYScpCm0xX2JldGEgPC0gbTEkY29lZmZpY2llbnRzCmxpbmVzKHNlcV9hcmVhLCBtMV9iZXRhWzFdICsgbTFfYmV0YVsyXSAqIHNlcV9hcmVhLCBjb2wgPSAncmVkJywgbHdkID0gMikKCnBsb3RfcmVudHNxbV9hcmVhKCdmKGFyZWEpID0gMS9hcmVhJykKbTJfYmV0YSA8LSBtMiRjb2VmZmljaWVudHMKbGluZXMoc2VxX2FyZWEsIG0yX2JldGFbMV0gKyBtMl9iZXRhWzJdICogMSAvIHNlcV9hcmVhLCBjb2wgPSAncmVkJywgbHdkID0gMikKCnBsb3RfcmVzaWR1YWxzIDwtIGZ1bmN0aW9uKG1vZGVsKSB7CiAgcGxvdCgKICAgIG11bmljaF9yZW50X2luZGV4JGFyZWEsCiAgICBtb2RlbCRyZXNpZHVhbHMsCiAgICBtYWluID0gJ3Jlc2lkdWFscycsCiAgICB4bGFiID0gJ2FyZWEgaW4gc3FtJywKICAgIHlsYWIgPSAncmVzaWR1YWxzJwogICkKICBhYmxpbmUoMCwgMCwgY29sID0gJ3JlZCcsIGx3ZCA9IDIpCn0KCnBsb3RfcmVzaWR1YWxzKG0xKQpwbG90X3Jlc2lkdWFscyhtMikKCnBsb3RfYXZfcmVzaWR1YWxzIDwtIGZ1bmN0aW9uKG1vZGVsKSB7CiAgYXZfcmVzaWR1YWxzIDwtIHRhcHBseSgKICAgIG1vZGVsJHJlc2lkdWFscywgbXVuaWNoX3JlbnRfaW5kZXgkYXJlYSwKICAgIG1lYW4KICApCiAgcGxvdCgKICAgIHVuaXF1ZShtdW5pY2hfcmVudF9pbmRleCRhcmVhKSwKICAgIGF2X3Jlc2lkdWFscywKICAgIG1haW4gPSAnYXZlcmFnZSByZXNpZHVhbHMnLAogICAgeGxhYiA9ICdhcmVhIGluIHNxbScsCiAgICB5bGFiID0gJ3Jlc2lkdWFscycKICApCiAgYWJsaW5lKDAsIDAsIGNvbCA9ICdyZWQnLCBsd2QgPSAyKQp9CgpwbG90X2F2X3Jlc2lkdWFscyhtMSkKcGxvdF9hdl9yZXNpZHVhbHMobTIpCmBgYAoKCiMjIEV4YW1wbGUgMy40IE11bmljaCBSZW50IEluZGV4IOKAlCBQb2x5bm9taWFsIFJlZ3Jlc3Npb24KCkZvciBwb2x5bm9taWFsIHJlZ3Jlc3Npb25zIHdlIGFkanVzdCB0d28gZGlmZmVyZW50IG1vZGVsLiBPZiBzZWNvbmQtb3JkZXIgKCoqTW9kZWwgMyAoJE0zJCkqKik6CgokJApcbWF0aGJie0V9KFx0ZXh0e3JlbnRxc219KSA9IFxiZXRhX3swfSArIFxiZXRhX3sxfSBcdGV4dHthcmVhfSArIFxiZXRhX3syfSBcdGV4dHthcmVhfV57Mn0KJCQKCkFuZCB0aGlyZC1vcmRlciAoKipNb2RlbCA0ICgkTTQkKSoqKToKCiQkClxtYXRoYmJ7RX0oXHRleHR7cmVudHFzbX0pID0gXGJldGFfezB9ICsgXGJldGFfezF9IFx0ZXh0e2FyZWF9ICsgXGJldGFfezJ9IFx0ZXh0e2FyZWF9XnsyfSArIFxiZXRhX3szfSBcdGV4dHthcmVhfV57M30KJCQKCmBgYHtyLCBmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9MTIsIGZpZy5hbGlnbj0nY2VudGVyJ30KcGFyKG1mcm93ID0gYygyLCAyKSkKCm0zIDwtIGxtKHJlbnRzcW0gfiBhcmVhICsgSShhcmVhXjIpLCBkYXRhID0gbXVuaWNoX3JlbnRfaW5kZXgpCnN1bW1hcnkobTMpCgptNCA8LSBsbShyZW50c3FtIH4gYXJlYSArIEkoYXJlYV4yKSArIEkoYXJlYV4zKSwgZGF0YSA9IG11bmljaF9yZW50X2luZGV4KQpzdW1tYXJ5KG00KQoKcGxvdF9yZW50c3FtX2FyZWEoJ3JlbnQgcGVyIHNxbSB2cy4gYXJlYScpCm0zX2JldGEgPC0gbTMkY29lZmZpY2llbnRzCmxpbmVzKHNlcV9hcmVhLCBtM19iZXRhWzFdICsgbTNfYmV0YVsyXSAqIHNlcV9hcmVhICsgbTNfYmV0YVszXSAqIHNlcV9hcmVhIF4gMiwgY29sID0gJ3JlZCcsIGx3ZCA9IDIpCgpwbG90X3JlbnRzcW1fYXJlYSgnZihhcmVhKSA9IDEvYXJlYScpCm00X2JldGEgPC0gbTQkY29lZmZpY2llbnRzCmxpbmVzKHNlcV9hcmVhLCBtNF9iZXRhWzFdICsgbTRfYmV0YVsyXSAqIHNlcV9hcmVhICsgbTRfYmV0YVszXSAqIHNlcV9hcmVhIF4gMiArIG00X2JldGFbNF0gKiBzZXFfYXJlYSBeIDMsIGNvbCA9ICdyZWQnLCBsd2QgPSAyKQoKcGxvdF9yZXNpZHVhbHMobTMpCnBsb3RfcmVzaWR1YWxzKG00KQpgYGAKCiMjIEV4YW1wbGUgMy41IE11bmljaCBSZW50IEluZGV4IOKAlCBBZGRpdGl2ZSBNb2RlbHMKCldlIGRlZmluZSB0aGUgZm9sbG93aW5nIGFkaXRpdmUgbW9kZWxzLiBBZGRpdGl2ZSBtb2RlbCAxOgoKJCQKXG1hdGhjYWx7RX0oXHRleHR7cmVudHNxbX1fe2l9KSA9IFxiZXRhX3swfSArIFxiZXRhX3sxfSBcY2RvdCBcZnJhY3sxfXtcdGV4dHthcmVhfV97aX19ICsgXGJldGFfezJ9IFxjZG90IFx0ZXh0e3llYXJjfV97aX0KJCQKCkFuZCBhIHNlY29uZCBvbmUgd2l0aCAkXHRleHR7eWVhcmN9JCBhcyBhIHBvbGlub21pYWwgb2YgZGVncmVlIDM6CgokJApcbWF0aGNhbHtFfShcdGV4dHtyZW50c3FtfV97aX0pID0gXGJldGFfezB9ICsgXGJldGFfezF9IFxjZG90IFxmcmFjezF9e1x0ZXh0e2FyZWF9X3tpfX0gKyBcYmV0YV97Mn0gXGNkb3QgXHRleHR7eWVhcmN9X3tpfSArIFxiZXRhX3szfSBcY2RvdCB7XHRleHR7eWVhcmN9X3tpfX1eezJ9ICsgXGJldGFfezR9IFxjZG90IHtcdGV4dHt5ZWFyY31fe2l9fV57M30KJCQKCmBgYHtyfQphZGRpdGl2ZV9tMSA8LSBsbShyZW50c3FtIH4gSSgxL2FyZWEpICsgeWVhcmMsIGRhdGEgPSBtdW5pY2hfcmVudF9pbmRleCkKc3VtbWFyeShhZGRpdGl2ZV9tMSkKCmFkZGl0aXZlX20yIDwtIGxtKHJlbnRzcW0gfiBJKDEvYXJlYSkgKyB5ZWFyYyArIEkoeWVhcmMgXiAyKSsgSSh5ZWFyYyBeIDMpLCBkYXRhID0gbXVuaWNoX3JlbnRfaW5kZXgpCnN1bW1hcnkoYWRkaXRpdmVfbTIpCmBgYAoKQSB0aGlyZCBtb2RlbCBpcyBzcGVjaWZpZWQgdXNpbmcgdGhlIHRydW5jYXRlZCB5ZWFyIG9mIGNvbnN0cnVjdGlvbiAoJFx0ZXh0e3llYXJjfSAtIDE5MDAkKToKCmBgYHtyfQptdW5pY2hfcmVudF9pbmRleCR5ZWFyY28gPC0gbXVuaWNoX3JlbnRfaW5kZXgkeWVhcmMgLSAxOTAwCmFkZGl0aXZlX20zIDwtIGxtKHJlbnRzcW0gfiBJKDEvYXJlYSkgKyB5ZWFyY28gKyBJKHllYXJjbyBeIDIpKyBJKHllYXJjbyBeIDMpLCBkYXRhID0gbXVuaWNoX3JlbnRfaW5kZXgpCnN1bW1hcnkoYWRkaXRpdmVfbTMpCmBgYAoKQXMgaW4gdGhlIGJvb2ssIHdlIHNob3cgdGhlIGVmZmVjdCBvZiB5ZWFyIG9mIGNvbnN0cnVjdGlvbiBpbiB0aGUgcG9seW5vbWlhbCBvZiBkZWdyZWUgMzoKCmBgYHtyLCBmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9NCwgZmlnLmFsaWduPSdjZW50ZXInfQpwYXIobWZyb3cgPSBjKDEsIDIpKQpiZXRhX2FkZF9tMiA8LSBhZGRpdGl2ZV9tMiRjb2VmZmljaWVudHMKCnllYXJjX3NlcSA8LSBzZXEobWluKG11bmljaF9yZW50X2luZGV4JHllYXJjKSwgbWF4KG11bmljaF9yZW50X2luZGV4JHllYXJjKSwgbGVuZ3RoLm91dCA9IDEwMCkKCnBsb3RfZWZmZWN0IDwtIGZ1bmN0aW9uKGJldGEsIHllYXJfc2VxLCB0aXRsZSl7CiAgcGxvdCgKICAgIHllYXJfc2VxLAogICAgYmV0YVszXSAqIHllYXJfc2VxICsgYmV0YVs0XSAqIHllYXJfc2VxXjIgKyBiZXRhWzVdICogeWVhcl9zZXFeMywKICAgIG1haW4gPSB0aXRsZSwKICAgIHhsYWIgPSAneWVhciBvZiBjb25zdHJ1Y3Rpb24nLAogICAgeWxhYiA9ICdlZmZlY3Qgb2YgeWVhciBvZiBjb25zdHJ1Y3Rpb24nLAogICAgdHlwZSA9ICdsJywKICAgIGNvbCA9ICdkYXJrYmx1ZScsCiAgICBsd2QgPSAyCiAgKQp9CgpwbG90X2VmZmVjdChiZXRhX2FkZF9tMiwgeWVhcmNfc2VxLCAnZWZmZWN0IG9mIHllYXIgb2YgY29uc3RydWN0aW9uJykKCmJldGFfYWRkX20zIDwtIGFkZGl0aXZlX20zJGNvZWZmaWNpZW50cwp0cnVuY2F0ZV95ZWFyY19zZXEgPC0gc2VxKG1pbihtdW5pY2hfcmVudF9pbmRleCR5ZWFyY28pLCBtYXgobXVuaWNoX3JlbnRfaW5kZXgkeWVhcmNvKSwgbGVuZ3RoLm91dCA9IDEwMCkKcGxvdF9lZmZlY3QoYmV0YV9hZGRfbTMsIHRydW5jYXRlX3llYXJjX3NlcSwgJ2VmZmVjdCBvZiB5ZWFyIG9mIGNvbnN0cnVjdGlvbiwgY29kZWQgZnJvbSAxOCB0byA5NicpCmBgYAoKQSBuZXcgbW9kZWwgaXMgZGVmaW5lZCB1c2luZyB0aGUgb3J0aG9nb25hbCBwb2x5bm9taWFsIG9mIHllYXIgb2YgY29uc3RydWN0aW9uOgoKYGBge3J9Cm11bmljaF9yZW50X2luZGV4JGFyZWFpbnYgPC0gKDEgLyBtdW5pY2hfcmVudF9pbmRleCRhcmVhKSAtIG1lYW4oMSAvIG11bmljaF9yZW50X2luZGV4JGFyZWEpCnBvbHlfeWVhciA8LSBwb2x5KG11bmljaF9yZW50X2luZGV4JHllYXJjbyAtIG1lYW4obXVuaWNoX3JlbnRfaW5kZXgkeWVhcmNvKSwgMykKbXVuaWNoX3JlbnRfaW5kZXgkeWVhcmNjIDwtIHBvbHlfeWVhclssIDFdCm11bmljaF9yZW50X2luZGV4JHllYXJjYzIgPC0gcG9seV95ZWFyWywgMl0KbXVuaWNoX3JlbnRfaW5kZXgkeWVhcmNjMyA8LSBwb2x5X3llYXJbLCAzXQphZGRpdGl2ZV9tNCA8LSBsbShyZW50c3FtIH4gYXJlYWludiArIHllYXJjYyArIHllYXJjYzIgKyB5ZWFyY2MzLCBkYXRhID0gbXVuaWNoX3JlbnRfaW5kZXgpCnN1bW1hcnkoYWRkaXRpdmVfbTQpCmBgYAoKTm93IHdlIGNhbGN1bGF0ZSB0aGUgcGFydGlhbCByZXNpZHVhbHMgZm9yICRcdGV4dHthcmVhfSQgYW5kICRcdGV4dHt5ZWFyY30kIGRlZmluZSBieToKCiQkClxoYXR7XGVwc2lsb259X3tcdGV4dHthcmVhfSwgaX0gPSBcdGV4dHtyZW50c3FtfV97aX0gLSBcaGF0e1xiZXRhX3swfX0gLSBcaGF0e2Z9X3syfShcdGV4dHt5ZWFyY31fe2l9KQokJAoKV2l0aCB0aGUgZWZmZWN0ICRmKFx0ZXh0e3llYXJjfV97aX0pJDoKCiQkClxoYXR7Zn1fezJ9KFx0ZXh0e3llYXJjfSkgPSBcYmV0YV97Mn0gXGNkb3QgXHRleHR7eWVhcmN9X3tpfSArIFxiZXRhX3szfSBcY2RvdCB7XHRleHR7eWVhcmN9X3tpfX1eezJ9ICsgXGJldGFfezR9IFxjZG90IHtcdGV4dHt5ZWFyY31fe2l9fV57M30KJCQKCkFuZDoKCiQkClxoYXR7XGVwc2lsb259X3tcdGV4dHt5ZWFyY30sIGl9ID0gXHRleHR7cmVudHNxbX1fe2l9IC0gXGhhdHtcYmV0YV97MH19IC0gXGhhdHtmfV97MX0oXHRleHR7YXJlYX1fe2l9KQokJAoKV2l0aCB0aGUgZWZmZWN0ICRmKFx0ZXh0e3llYXJjfV97aX0pJDoKCiQkClxoYXR7Zn1fezF9KFx0ZXh0e3llYXJjfSkgPSBcYmV0YV97MX0gXGNkb3QgXGZyYWN7MX17XHRleHR7YXJlYX1fe2l9fQokJAoKYGBge3IsIGZpZy53aWR0aD0xMCwgZmlnLmhlaWdodD00LCBmaWcuYWxpZ249J2NlbnRlcid9CnBhcihtZnJvdyA9IGMoMSwgMikpCgpiZXRhX2FkZF9tNCA8LSBhZGRpdGl2ZV9tNCRjb2VmZmljaWVudHMKCmFyZWFfc2VxIDwtIHNlcShtaW4obXVuaWNoX3JlbnRfaW5kZXgkYXJlYSksIG1heChtdW5pY2hfcmVudF9pbmRleCRhcmVhKSwgbGVuZ3RoLm91dCA9IDEwMCkKCmFyZWFfZWZmZWN0IDwtIGZ1bmN0aW9uKGFyZWEpewogIGludl9hcmVhIDwtIDEgLyBhcmVhCiAgY2VudF9hcmVhIDwtIGludl9hcmVhIC0gbWVhbihpbnZfYXJlYSkKICBiZXRhX2FkZF9tNFsyXSAqIGNlbnRfYXJlYQp9Cgp5ZWFyY19lZmZlY3QgPC0gZnVuY3Rpb24oeWVhcmMpewogIHllYXJjX2NlbnRlciA8LSB5ZWFyYyAtIG1lYW4oeWVhcmMpCiAgcG9seV95ZWFyYyA8LSBwb2x5KHllYXJjX2NlbnRlciwgMykKICBiZXRhX2FkZF9tNFszXSAqIHBvbHlfeWVhcmNbLCAxXSArIGJldGFfYWRkX200WzRdICogcG9seV95ZWFyY1ssIDJdICsgYmV0YV9hZGRfbTRbNV0gKiBwb2x5X3llYXJjWywgM10KfQoKcmVzaWR1YWxfYXJlYSA8LSBtdW5pY2hfcmVudF9pbmRleCRyZW50c3FtIC0gYmV0YV9hZGRfbTRbMV0gLSB5ZWFyY19lZmZlY3QobXVuaWNoX3JlbnRfaW5kZXgkeWVhcmMpCnBsb3QoCiAgbXVuaWNoX3JlbnRfaW5kZXgkYXJlYSwKICByZXNpZHVhbF9hcmVhLAogIG1haW4gPSAnZWZmZWN0IG9mIGFyZWEnLAogIHhsYWIgPSAnYXJlYSBpbiBzcW0nLAogIHlsYWIgPSAnZWZmZWN0IG9mIGFyZWEgLyBwYXJ0aWFsIHJlc2lkdWFscycKKQpsaW5lcyhhcmVhX3NlcSwgYXJlYV9lZmZlY3QoYXJlYV9zZXEpLCBjb2wgPSAncmVkJywgbHdkID0gMikKCnJlc2lkdWFsX3llYXIgPC0gbXVuaWNoX3JlbnRfaW5kZXgkcmVudHNxbSAtIGJldGFfYWRkX200WzFdIC0gYXJlYV9lZmZlY3QobXVuaWNoX3JlbnRfaW5kZXgkYXJlYSkKcGxvdCgKICBtdW5pY2hfcmVudF9pbmRleCR5ZWFyYywKICByZXNpZHVhbF95ZWFyLAogIG1haW4gPSAnZWZmZWN0IG9mIHllYXIgb2YgY29uc3RydWN0aW9uJywKICB4bGFiID0gJ3llYXIgb2YgY29uc3RydWN0aW9uJywKICB5bGFiID0gJ2VmZmVjdCBvZiB5ZWFyIG9mIGNvbnN0cnVjdGlvbiAvIHBhcnRpYWwgcmVzaWR1YWxzJwopCmxpbmVzKHllYXJjX3NlcSwgeWVhcmNfZWZmZWN0KHllYXJjX3NlcSksIGNvbCA9ICdyZWQnLCBsd2QgPSAyKQpgYGAKCiMjIEV4YW1wbGUgMy42IE11bmljaCBSZW50IEluZGV4IOKAlCBFZmZlY3QgQ29kaW5nCgpXZSBhZGp1c3QgdGhlIHJlZ3Jlc3Npb24gbW9kZWwgd2l0aCB0aGUgY29kaW5nIHZhbHVlcyBvZiBsb2NhdGlvbjoKCiQkClxtYXRoY2N7RX0oXHRleHR7cmVudHNxbX1fe2l9KSA9IFxiZXRhX3swfSArIFxiZXRhX3sxfSBcY2RvdCBcdGV4dHtnbG9jYXRpb259X3sxfSArIFxiZXRhX3syfSBcY2RvdCBcdGV4dHt0bG9jYXRpb259X3sxfQokJAoKYGBge3J9Cm11bmljaF9yZW50X2luZGV4JGdsb2NhdGlvbiA8LSAwCm11bmljaF9yZW50X2luZGV4JGdsb2NhdGlvblttdW5pY2hfcmVudF9pbmRleCRsb2NhdGlvbiA9PSAyXSA8LSAxCm11bmljaF9yZW50X2luZGV4JGdsb2NhdGlvblttdW5pY2hfcmVudF9pbmRleCRsb2NhdGlvbiA9PSAxXSA8LSAtMQoKbXVuaWNoX3JlbnRfaW5kZXgkdGxvY2F0aW9uIDwtIDAKbXVuaWNoX3JlbnRfaW5kZXgkdGxvY2F0aW9uW211bmljaF9yZW50X2luZGV4JGxvY2F0aW9uID09IDNdIDwtIDEKbXVuaWNoX3JlbnRfaW5kZXgkdGxvY2F0aW9uW211bmljaF9yZW50X2luZGV4JGxvY2F0aW9uID09IDFdIDwtIC0xCgpjb2RfbW9kZWwgPC0gbG0oCiAgcmVudHNxbSB+IGdsb2NhdGlvbiArIHRsb2NhdGlvbiwKICBkYXRhID0gbXVuaWNoX3JlbnRfaW5kZXgKKQpzdW1tYXJ5KGNvZF9tb2RlbCkKYGBgCgojIyBFeGFtcGxlIDMuNyBNdW5pY2ggUmVudCBJbmRleCDigJQgSW50ZXJhY3Rpb24gd2l0aCBRdWFsaXR5IG9mIEtpdGNoZW4KCldlIGltcG9ydCB0aGUgbmV3IHVwZGF0ZWQgbW9kZWwgb2YgMjAwMSwgYXZhaWxhYmxlIGF0IFtvZmZpY2lhbCBwYWdlIG9mIHRoZSBkYXRhc2V0XShodHRwczovL3d3dy51bmktZ29ldHRpbmdlbi5kZS9lbi81NTE2MjUuaHRtbCkgYXMgYWxsIHRoZSBkYXRhc2V0czoKCmBgYHtyfQptdW5pY2hfcmVudF91cmxfMDEgPC0gImh0dHBzOi8vd3d3LnVuaS1nb2V0dGluZ2VuLmRlL2RlL2RvY3VtZW50L2Rvd25sb2FkL2UwMGQwMzllOWMxZjI4YWI0MDRiMjdiMGIxNDkzMWYxLnJhdy9yZW50OThfMDAucmF3IgoKbXVuaWNoX3JlbnRfaW5kZXhfMDEgPC0gcmVhZC50YWJsZSgKICB1cmwobXVuaWNoX3JlbnRfdXJsXzAxKSwKICBoZWFkZXIgPSAxLAogIGNvbENsYXNzZXMgPSBjKAogICAgIm51bWVyaWMiLCAibnVtZXJpYyIsICJudW1lcmljIiwKICAgICJudW1lcmljIiwgImludGVnZXIiLCAiaW50ZWdlciIsCiAgICAiaW50ZWdlciIsICJpbnRlZ2VyIiwgImludGVnZXIiLAogICAgImludGVnZXIiLCAibnVtZXJpYyIsICJudW1lcmljIiwKICAgICJudW1lcmljIgogICkKKQoKc3VtbWFyeShtdW5pY2hfcmVudF9pbmRleF8wMSkKYGBgCgpBbmQgd2UgYWRqdXN0IHRoZSBtb2RlbDoKCiQkClxiZWdpbnttYXRyaXh9CiAgXHdpZGVoYXR7XHRleHR7cmVudHNxbX19X3tpfSA9ICYgXGJldGFfezB9ICsgXGJldGFfezF9IFx0ZXh0e2FyZWFpbnZjfV97aX0gKyBcYmV0YV97Mn0gXHRleHR7eWVhcmNvfV97aX0gKyBcYmV0YV97M30gXHRleHR7eWVhcmNvMn1fe2l9ICsgXGJldGFfezR9IFx0ZXh0e3llYXJjbzN9X3tpfSBcXAogICAmICsgXGJldGFfezV9IFx0ZXh0e3llYXIwMX1fe2l9ICsgXGJldGFfezZ9IFx0ZXh0e25raXRjaGVufV97aX0gKyBcYmV0YV97N30gXHRleHR7cGtpdGNoZW59X3tpfSBcXAogICAmICsgXGJldGFfezh9IFx0ZXh0e25raXRjaGVufV97aX0gXGNkb3QgXHRleHR7eWVhcjAxfV97aX0gKyBcYmV0YV97OX0gXHRleHR7cGtpdGNoZW59X3tpfSBcY2RvdCBcdGV4dHt5ZWFyMDF9X3tpfQpcZW5ke21hdHJpeH0KJCQKCmBgYHtyfQptdW5pY2hfcmVudF9pbmRleF8wMSRhcmVhaW52YyA8LSBtdW5pY2hfcmVudF9pbmRleF8wMSRpbnZhcmVhIC0gbWVhbihtdW5pY2hfcmVudF9pbmRleF8wMSRpbnZhcmVhKQptdW5pY2hfcmVudF9pbmRleF8wMSR5ZWFyY28gPC0gbXVuaWNoX3JlbnRfaW5kZXhfMDEkeWVhcmMgLSBtZWFuKG11bmljaF9yZW50X2luZGV4XzAxJHllYXJjKQpwb2x5X3llYXJfMDEgPC0gcG9seShtdW5pY2hfcmVudF9pbmRleF8wMSR5ZWFyY28sIDMpCm11bmljaF9yZW50X2luZGV4XzAxJHllYXJjbyA8LSBwb2x5X3llYXJfMDFbLCAxXQptdW5pY2hfcmVudF9pbmRleF8wMSR5ZWFyY28yIDwtIHBvbHlfeWVhcl8wMVssIDJdCm11bmljaF9yZW50X2luZGV4XzAxJHllYXJjbzMgPC0gcG9seV95ZWFyXzAxWywgM10KCmludGVyYWN0aW9uX21vZGVsIDwtIGxtKAogIHJlbnRzcW1fZXVybyB+IGFyZWFpbnZjICsgeWVhcmNvICsgeWVhcmNvMiArIHllYXJjbzMgKyB5ZWFyMDEgKyBua2l0Y2hlbiArIHBraXRjaGVuICsgbmtpdGNoZW4gKiB5ZWFyMDEgKyBwa2l0Y2hlbiAqIHllYXIwMSwKICBkYXRhID0gbXVuaWNoX3JlbnRfaW5kZXhfMDEKKQoKc3VtbWFyeShpbnRlcmFjdGlvbl9tb2RlbCkKYGBgCgojIyBFeGFtcGxlIDMuOCBNdW5pY2ggUmVudCBJbmRleCDigJQgSW50ZXJhY3Rpb24gQmV0d2VlbiBMaXZpbmcgQXJlYSBhbmQgTG9jYXRpb24KClRoZSBtb2RlbCB0byBhZGp1c3QgaXVzIGRlZmluZWQgYXM6CgokJApcbWF0aGJie0V9KFx0ZXh0e3JlbnRzcW19X3tpfSkgPSBcYmV0YV97MH0gKyBcYmV0YV97MX0gXGNkb3QgXHRleHR7dGxvY2F0aW9ufSArIFxiZXRhX3syfSBcY2RvdCBcb3ZlcmxpbmV7XGZyYWN7MX17XHRleHR7YXJlYX19fSArIFxiZXRhezN9IFxjZG90IFx0ZXh0e2FyZWF9IFxjZG90IFx0ZXh0e3Rsb2NhdGlvbn0KJCQKCldlIHVzZSB0aGUgZGF0YXNldCBmb3IgYXZlcmFnZSBvciB0b3AgbG9jYXRpb24gb25seSBhbmQgY2hhbmdlIHRoZSBkdW1teSB2YXJpYWJsZSAkXHRleHR7dGxvY2F0aW9ufSQgdG8gJDAkIChhdmVyYWdlKSBhbmQgJDEkICh0b3ApLiBGaW5hbGx5LCB3ZSB2aXN1YWxpemUgdGhlIGVmZmVjdCBhcyBkZXNjcmliZWQgaW4gdGhlIGJvb2suCgpgYGB7ciwgZmlnLndpZHRoPTEwLCBmaWcuaGVpZ2h0PTQsIGZpZy5hbGlnbj0nY2VudGVyJ30KcGFyKG1mcm93ID0gYygxLCAyKSkKYXRfcmVudCA8LSBzdWJzZXQobXVuaWNoX3JlbnRfaW5kZXgsIGxvY2F0aW9uID09IDEgfCBsb2NhdGlvbiA9PSAzKQphdF9yZW50JHRsb2NhdGlvblthdF9yZW50JGxvY2F0aW9uID09IDFdIDwtIDAKc3VtbWFyeShhdF9yZW50KQoKYXJlYV9sb2NhdGlvbl9tb2RlbCA8LSBsbSgKICByZW50c3FtIH4gdGxvY2F0aW9uICsgYXJlYWludiArIGFyZWE6dGxvY2F0aW9uLAogIGRhdGEgPSBhdF9yZW50CikKc3VtbWFyeShhcmVhX2xvY2F0aW9uX21vZGVsKQoKYmV0YV9hbCA8LSBhcmVhX2xvY2F0aW9uX21vZGVsJGNvZWZmaWNpZW50cwoKZjFfYXJlYSA8LSBmdW5jdGlvbihhcmVhKXsKICBhcmVhaW52YyA8LSAoMSAvIGFyZWEpIC0gbWVhbigxIC8gYXJlYSkKICBiZXRhX2FsWzNdICogYXJlYWludmMKfQoKZjJfYXJlYSA8LSBmdW5jdGlvbihhcmVhKXsKICBiZXRhX2FsWzRdICogYXJlYQp9CgpzbWFsbF9lZmZlY3QgPC0gZjFfYXJlYShhcmVhX3NlcSkKYmlnX2VmZmVjdCA8LSBiZXRhX2FsWzJdICsgZjFfYXJlYShhcmVhX3NlcSkgKyBmMl9hcmVhKGFyZWFfc2VxKQoKeWxpbSA8LSBjKG1pbihjKHNtYWxsX2VmZmVjdCwgYmlnX2VmZmVjdCkpLCBtYXgoYyhzbWFsbF9lZmZlY3QsIGJpZ19lZmZlY3QpKSkKCnBsb3QoCiAgYXJlYV9zZXEsCiAgc21hbGxfZWZmZWN0LAogIG1haW4gPSAnYXJlYSBlZmZlY3RzIGJhc2VkIG9uIGEgbGluZWFyIGludGVyYWN0aW9uJywKICB4bGFiID0gJ2FyZWEgaW4gc3FtJywKICB5bGFiID0gJ2FyZWEgZWZmZWN0cycsCiAgeWxpbSA9IHlsaW0sCiAgdHlwZSA9ICdsJywKICBjb2wgPSAnZGFya2JsdWUnLAogIGx3ZCA9IDIKKQpsaW5lcyhhcmVhX3NlcSwgYmlnX2VmZmVjdCwgY29sID0gJ3JlZCcsIGx3ZCA9IDIpCgpwbG90KAogIGFyZWFfc2VxLAogIGJldGFfYWxbMl0gKyBmMl9hcmVhKGFyZWFfc2VxKSwKICBtYWluID0gJ3ZhcnlpbmcgZWZmZWN0IG9mIHRvcCBsb2NhdGlvbicsCiAgeGxhYiA9ICdhcmVhIGluIHNxbScsCiAgeWxhYiA9ICdWYXJ5aW5nIGVmZmVjdCBvZiB0b3AgbG9jYXRpb24nLAogIHR5cGUgPSAnbCcsCiAgY29sID0gJ2RhcmtibHVlJywKICBsd2QgPSAyCikKYGBgCgojIyBFeGFtcGxlIDMuOSBPcnRob2dvbmFsIERlc2lnbgoKSnVzdCBmb3IgZXhlbXBsaWZ5LCB3ZSBkZWZpbmUgdGhlIG9ydGhvZ29uYWwgcHJvY2VzcyBkZXNjcmliZSBhcyBhbiBSIGZ1bmN0aW9uOgoKJCQKXHRpbGRle3h9XntqfSA9IHhee2p9IC0gXHdpZGV0aWxkZXtYfV97an0gXGxlZnQgKCB7XHdpZGV0aWxkZXtYfV97an19XntUfSBcd2lkZXRpbGRle1h9X3tqfSBccmlnaHQgKV57LTF9IHtcd2lkZXRpbGRle1h9X3tqfX1ee1R9IHhee2p9CiQkCgpgYGB7cn0Kb3J0aG9nb25hbF9kZXNpZ24gPC0gZnVuY3Rpb24oeCl7CiAgeF9vdXQgPC0gY2JpbmQoeFssIDFdIC0gbWVhbih4WywgMV0pLCBtYXRyaXgoMCwgbnJvdyA9IG5yb3coeCksIG5jb2wgPSBuY29sKHgpIC0gMSkpCiAgZm9yIChpIGluIDI6bmNvbCh4KSl7CiAgICB4X2hhdCA8LSBhcy5tYXRyaXgoeF9vdXRbLCAxOihpLTEpXSkKICAgIHhfaW52IDwtIHNvbHZlKHQoeF9oYXQpICUqJSB4X2hhdCkKICAgIHhfb3V0WywgaV0gPC0geFssIGldIC0geF9oYXQgJSolIHhfaW52ICUqJSB0KHhfaGF0KSAlKiUgeFssIGldCiAgfQogIHhfb3V0Cn0KYGBgCldlIHRlc3QgdGhpcyBmdW5jdGlvbiBvbiB0aGUgcG9seW5vbWlhbCBvZiB0aGUgJFx0ZXh0e3llYXJjfSQgdmFyaWFibGU6CgpgYGB7cn0KcG9seV90ZXN0IDwtIHBvbHkobXVuaWNoX3JlbnRfaW5kZXgkeWVhcmMsIDMsIHJhdyA9IFRSVUUpCm9ydGhvZ29uYWxfdGVzdCA8LSBvcnRob2dvbmFsX2Rlc2lnbihwb2x5X3Rlc3QpCgojIEFuZCBjaGVjayBpdApjKAogIG9ydGhvZ29uYWxfdGVzdFssIDFdICUqJSBvcnRob2dvbmFsX3Rlc3RbLCAyXSwKICBvcnRob2dvbmFsX3Rlc3RbLCAyXSAlKiUgb3J0aG9nb25hbF90ZXN0WywgM10sCiAgb3J0aG9nb25hbF90ZXN0WywgMV0gJSolIG9ydGhvZ29uYWxfdGVzdFssIDNdCikKYGBgCgpXaHkgZGlkIGl0IG5vdCByZXN1bHQ/CgojIyBFeGFtcGxlIDMuMTAgTXVuaWNoIFJlbnQgSW5kZXgg4oCUIENvbXBhcmlzb24gb2YgTW9kZWxzIFVzaW5nICR7Un1eezJ9JAoKV2UgY29tcGFyZSB0aGUgbW9kZWxzOiAkTTEkLCAkTTIkLCAkTTMkLCBhbmQgJE00JCB1c2luZyBpdHMgY29lZmZpY2llbnQgb2YgZGV0ZXJtaW5hdGlvbiAkUl57Mn0kCgpgYGB7cn0KY29lZmZfZGV0ZXIgPC0gZGF0YS5mcmFtZSgKICBgUiBTcXVhcmVkYCA9IGMoc3VtbWFyeShtMSkkci5zcXVhcmVkLCBzdW1tYXJ5KG0yKSRyLnNxdWFyZWQsIHN1bW1hcnkobTMpJHIuc3F1YXJlZCwgc3VtbWFyeShtNCkkci5zcXVhcmVkKSwKICByb3cubmFtZXMgPSBjKCJNMSIsICJNMiIsICJNMyIsICJNNCIpCikKY29lZmZfZGV0ZXIKYGBgCgojIyBFeGFtcGxlIDMuMTEgU2ltcGxlIExpbmVhciBSZWdyZXNzaW9uCgpGb3IgdGhlIG1vZGVsOgoKJCQKeV97aX0gPSBcYmV0YSB4X3tpfSArIFxlcHNpbG9uX3tpfQokJAoKV2UgY2FuIHZlcmlmeToKCiQkClxsZWZ0ICgge1hfe259fV57VH0gWF97bn0gXHJpZ2h0ICleey0xfSBcdG8gMAokJAoKQW5kOgoKJCQKXG1heF97aT0xLC4uLixufSB7eF97aX19XntUfSBcbGVmdCAoIHtYX3tufX1ee1R9IFhfe259IFxyaWdodCApXnstMX0geF97aX0gXHRvIDAgXHRleHR7IGZvciB9IG4gXHRvIFxpbmZ0eQokJAoKIyMjICR4X3tpfSA9IGkkOgoKJCQKXGxlZnQgKCB7WF97bn19XntUfSBYX3tufSBccmlnaHQgKV57LTF9ID0gXGxlZnQgKCAoMSwgMiwgMywgXGNkb3RzLCBuKSBcbGVmdCAoIFxiZWdpbnttYXRyaXh9CjEgXFwKMiBcXApcdmRvdHMgIFxcCm4KXGVuZHttYXRyaXh9IFxyaWdodCApIFxyaWdodCApXnstMX0gPSBcbGVmdCAoIDEgKyA0ICsgXGNkb3RzICsgbl4yIFxyaWdodCApXnstMX0gXHRvIDAKJCQKCkFuZDoKCiQkClxtYXhfe2k9MSwuLi4sbn0gKCAxLCAyLCBcY2RvdHMsIG4gKSBcbGVmdCAoIDEgKyA0ICsgXGNkb3RzICsgbl4yIFxyaWdodCApXnstMX0gXGxlZnQgKCBcYmVnaW57bWF0cml4fQoxIFxcCjIgXFwKXHZkb3RzICBcXApuClxlbmR7bWF0cml4fSBccmlnaHQgKSBcdG8gMCBcdGV4dHsgZm9yIH0gbiBcdG8gXGluZnR5CiQkCgojIyMgJHhfe2l9ID0gXGZyYWN7MX17aX0kOgoKJCQKXGxlZnQgKCB7WF97bn19XntUfSBYX3tufSBccmlnaHQgKV57LTF9ID0gXGxlZnQgKCBcbGVmdCAoMSwgMC41LCAwLjMzLCBcY2RvdHMsIFxmcmFjezF9e259IFxyaWdodCApIFxsZWZ0ICggXGJlZ2lue21hdHJpeH0KMSBcXAowLjUgXFwKXHZkb3RzICBcXApcZnJhY3sxfXtufQpcZW5ke21hdHJpeH0gXHJpZ2h0ICkgXHJpZ2h0ICleey0xfSA9IFxsZWZ0ICggXHN1bV97aT0xfV57bn0gXGZyYWN7MX17aV57Mn19IFxyaWdodCApXnstMX0gXG5vdHtcdG99IDAKJCQKCkFuZDoKCiQkClxtYXhfe2k9MSwuLi4sbn0gXGxlZnQgKCAxLCAwLjUsIFxjZG90cywgXGZyYWN7MX17bn0gXHJpZ2h0ICkgXGxlZnQgKCBcc3VtX3tpID0gMX1ee259IFxmcmFjezF9e2leezJ9fSBccmlnaHQgKV57LTF9IFxsZWZ0ICggXGJlZ2lue21hdHJpeH0KMSBcXAowLjUgXFwKXHZkb3RzICBcXApcZnJhY3sxfXtufQpcZW5ke21hdHJpeH0gXHJpZ2h0ICkgXG5vdHtcdG99IDAgXHRleHR7IGZvciB9IG4gXHRvIFxpbmZ0eQokJAoKIyMjICR4X3tpfSA9IFxmcmFjezF9e1xzcXJ0e2l9fSQ6CgokJApcbGVmdCAoIHtYX3tufX1ee1R9IFhfe259IFxyaWdodCApXnstMX0gPSBcbGVmdCAoIFxzdW1fe2kgPSAxfV57bn0gXGZyYWN7MX17aX0gXHJpZ2h0ICleey0xfSBcdG8gMAokJAoKQW5kOgoKJCQKXG1heF97aT0xLC4uLixufSBcbGVmdCAoIDEsIFxmcmFjezF9e1xzcXJ0ezJ9fSwgXGNkb3RzLCBcZnJhY3sxfXtcc3FydHtufX0gXHJpZ2h0ICkgXGxlZnQgKCBcc3VtX3tpID0gMX1ee259IFxmcmFjezF9e2l9IFxyaWdodCApXnstMX0gXGJlZ2lue21hdHJpeH0KMSBcXApcZnJhY3sxfXtcc3FydHsyfX0gXFwKXHZkb3RzIFxcClxmcmFjezF9e1xzcXJ0e259fQpcZW5ke21hdHJpeH0gXHRvIDAKJCQKCiMjIEV4YW1wbGUgMy4xMiBNdW5pY2ggUmVudCBJbmRleCDigJQgSHlwb3RoZXNpcyBUZXN0aW5nCgpUaGUgcmVncmVzc2lvbiBtb2RlbCB0byBhZGp1c3QgaXM6CgokJApcYmVnaW57bWF0cml4fQpcdGV4dHtyZW50c3FtfV97aX0gPSAmIFxiZXRhX3swfSArIFxiZXRhX3sxfSBcdGV4dHthcmVhaW52Y31fe2l9ICsgXGJldGFfezJ9IFx0ZXh0e3llYXJjb31fe2l9ICsgXGJldGFfezN9IFx0ZXh0e3llYXJjbzJ9X3tpfSArIFxiZXRhX3s0fSBcdGV4dHt5ZWFyY28zfV97aX0gXFwKJiArIFxiZXRhX3s1fSBcdGV4dHtua2l0Y2hlbn1fe2l9ICsgXGJldGFfezZ9IFx0ZXh0e3BraXRjaGVufV97aX0gKyBcYmV0YV97N30gXHRleHR7eWVhcjAxfV97aX0gKyBcZXBzaWxvbl97aX0KXGVuZHttYXRyaXh9CiQkCgpgYGB7cn0KaHlwb3RoZXNpc19tb2RlbCA8LSBsbSgKICByZW50c3FtX2V1cm8gfiBhcmVhaW52YyArIHllYXJjbyArIHllYXJjbzIgKyB5ZWFyY28zICsgbmtpdGNoZW4gKyBwa2l0Y2hlbiArIHllYXIwMSwKICBkYXRhID0gbXVuaWNoX3JlbnRfaW5kZXhfMDEKKQpzdW1tYXJ5KGh5cG90aGVzaXNfbW9kZWwpCmBgYAoKU28gd2Ugd2FudCB0byB0ZXN0IHRoZSBoeXBvdGhlc2lzOgoKJCQKXGJlZ2lue21hdHJpeH0KSF97MH0gOiBcYmV0YV97N30gPSAwICYgXHRleHR7YWdhaW5zdH0gJiBIX3sxfSA6IFxiZXRhX3s3fSBcbmVxIDAKXGVuZHttYXRyaXh9CiQkCgpBYm91dCB0aGUgc2lnbmlmaWNhbmNlIG9mIHRoZSBmb2xsb3ctdXAgbWVhc3VyZS4KCiQkClxiZWdpbnttYXRyaXh9CkhfezB9IDogXGxlZnQgKCBcYmVnaW57bWF0cml4fQpcYmV0YV97NX0gXFwKXGJldGFfezZ9ClxlbmR7bWF0cml4fSBccmlnaHQgKSA9IFxsZWZ0ICggXGJlZ2lue21hdHJpeH0KMCBcXAowClxlbmR7bWF0cml4fSBccmlnaHQgKSAmIFx0ZXh0e2FnYWluc3R9ICYgSF97MX0gOiBcbGVmdCAoIFxiZWdpbnttYXRyaXh9ClxiZXRhX3s1fSBcXApcYmV0YV97Nn0KXGVuZHttYXRyaXh9IFxyaWdodCApIFxuZXEgXGxlZnQgKCBcYmVnaW57bWF0cml4fQowIFxcCjAKXGVuZHttYXRyaXh9IFxyaWdodCApClxlbmR7bWF0cml4fQokJAoKQWJvdXQgdGhlIHNpZ25pZmljYW5jZSBvZiB0aGUga2l0Y2hlbiB2YXJpYWJsZS4gQW5kIGEgZmluYWwgdGVzdCBhYm91dCB0aGUgc2lnbmlmaWNhbmNlIG9mIGRpZmZlcmVudGlhdGUgYmV0d2VlbiB0aGVzZSB0d28gdHlwZSBvZiBraXRjaGVuczoKCiQkClxiZWdpbnttYXRyaXh9CkhfezB9IDogXGJldGFfezV9IC0gXGJldGFfezZ9ID0gMCAmIFx0ZXh0e2FnYWluc3R9ICYgSF97MX0gOiBcYmV0YV97NX0gLSBcYmV0YV97Nn0gXG5lcSAwClxlbmR7bWF0cml4fQokJAoKIyMgRXhhbXBsZSAzLjEzIE11bmljaCBSZW50IEluZGV4IOKAlCBTdGFuZGFyZCBPdXRwdXQgYW5kIEh5cG90aGVzaXMgVGVzdHMKCkZvciB0aGUgZmlyc3QgdGVzdCAoJEhfezB9IDogXGJldGFfezd9ID0gMCQpLCB3ZSBjYWxjdWxhdGU6CgokJAp0X3s3fSA9IFxmcmFje1xoYXR7XGJldGF9X3s3fX17e1x3aWRlaGF0e1x0ZXh0e1Zhcn0oXGhhdHtcYmV0YX1fezd9KX19XnsxLzJ9fSBcc2ltIHRfezEgLSBcYWxwaGEgLyAyfSAobiAtIHApCiQkCgpUaGUgY292YXJpYW5jZSBtYXRyaXggaXM6CgpgYGB7cn0KY292X21hdHJpeCA8LSB2Y292KGh5cG90aGVzaXNfbW9kZWwpCmNvdl9tYXRyaXgKYGBgCgpXZSBjaGVjayB0aGUgJFx0ZXh0e1Zhcn0oXGhhdHtcYmV0YX1fezd9KSQgaXMgdGhlIGxhc3QgdmFsdWUuCgpgYGB7cn0KdmFyXzcgPC0gY292X21hdHJpeFs4LCA4XQp0XzcgPC0gaHlwb3RoZXNpc19tb2RlbCRjb2VmZmljaWVudHNbOF0gLyBzcXJ0KHZhcl83KQp0XzcKYGBgCgpBbmQgdGhpcyBjb2luY2lkZXMgd2l0aCB0aGUgdmFsdWUgZ2l2ZW4gaW4gdGhlIGBzdW1tYXJ5YC4KCkZvciB0aGUgc2Vjb25kIHRlc3QsIHdlIGNvbXB1dGUgdGhlICRGJCB2YWx1ZSBhczoKCiQkCkYgPSBcZnJhY3sxfXtyfSB7XGhhdHtcYmV0YX19XntUfSBcd2lkZWhhdHtcdGV4dHtDb3Z9IFxsZWZ0ICggXGhhdHtcYmV0YX0gXHJpZ2h0ICl9XnstMX0gXGhhdHtcYmV0YX0gXHNpbSBGX3sxLCBuIC0gcH0KJCQKCldpdGggdGhlICRcaGF0e1xiZXRhfSQgZGVmaW5lIGluIHRoZSB0ZXN0ICRcbGVmdCAoIFxiZWdpbnttYXRyaXh9ClxoYXR7XGJldGFfezV9fSBcXApcaGF0e1xiZXRhX3s2fX0KXGVuZHttYXRyaXh9IFxyaWdodCApJCwgdGhlbiAkcj0yJCBhbmQgJG4gLSBwID0gbiAtIDgkOgoKYGBge3J9CmRmIDwtIG5yb3cobXVuaWNoX3JlbnRfaW5kZXhfMDEpIC0gOApyIDwtIDIKYmV0YV9oYXQgPC0gYXMubWF0cml4KGh5cG90aGVzaXNfbW9kZWwkY29lZmZpY2llbnRzWzY6N10pCmZfNTYgPC0gKDEvMikgKiB0KGJldGFfaGF0KSAlKiUgc29sdmUoY292X21hdHJpeFs2OjcsIDY6N10pICUqJSBiZXRhX2hhdApmXzU2CmBgYAoKV2UgZ2V0IGV4cGVjdGVkICRGJDoKCmBgYHtyfQpmX2NyaXQgPC0gcWYocCA9IDAuOTUsIHIsIGRmKQpmX2NyaXQKYGBgCgpBbmQgd2UgcmVqZWN0IHRoZSBudWxsIGh5cG90aGVzaXMuCgpUaGUgdGhpcmQgYW5kIGZpbmFsIHRlc3QgY29tcGFyaW5nICRcYmV0YV97NX0kIGFuZCAkXGJldGFfezZ9JCBiZXR3ZWVuIGVhY2ggb3RoZXIuIFdlIG5lZWQ6CgokJApcd2lkZWhhdHtcdGV4dHtWYXJ9KFxoYXR7XGJldGF9X3s1fSAtIFxoYXR7XGJldGF9X3s2fSl9ID0gXHdpZGVoYXR7XHRleHR7VmFyfShcaGF0e1xiZXRhfV97NX0pfSArIFx3aWRlaGF0e1x0ZXh0e1Zhcn0oXGhhdHtcYmV0YX1fezZ9KX0gLSAyIFxjZG90IFx3aWRlaGF0e1x0ZXh0e0Nvdn0oXGhhdHtcYmV0YX1fezV9LCBcaGF0e1xiZXRhfV97Nn0pfQokJAoKVXNpbmcgdGhlICRcd2lkZWhhdHtcdGV4dHtDb3Z9KFxoYXR7XGJldGF9KX0kIG1hdHJpeDoKCmBgYHtyfQpjb3ZfNV82IDwtIGNvdl9tYXRyaXhbNiwgN10KdmFyXzUgPC0gY292X21hdHJpeFs2LCA2XQp2YXJfNiA8LSBjb3ZfbWF0cml4WzcsIDddCnZhcl81XzYgPC0gdmFyXzUgKyB2YXJfNiAtIDIgKiBjb3ZfNV82CmZfNV82IDwtIChoeXBvdGhlc2lzX21vZGVsJGNvZWZmaWNpZW50c1s2XSAtIGh5cG90aGVzaXNfbW9kZWwkY29lZmZpY2llbnRzWzddKV4yIC8gdmFyXzVfNgpmXzVfNgpgYGAKCkFuZCB0aGlzICRGJCBjcml0aWNhbCwgdXNpbmcgJHI9MSQ6CgpgYGB7cn0KZl9jcml0IDwtIHFmKHAgPSAwLjk1LCAxLCBkZikKZl9jcml0CmBgYAoKU28sIGluIHRoaXMgY2FzZSwgd2UgZG8gbm90IHJlamVjdCB0aGUgbnVsbCBoeXBvdGhlc2lzLgoKIyMgRXhhbXBsZSAzLjE0IE11bmljaCBSZW50IEluZGV4IOKAlCBDb25maWRlbmNlIEludGVydmFscwoKVGhlIGNvbmZpZGVuY2UgaW50ZXJ2YWwgZm9yICRcYmV0YV97N30kIGlzIG9idGFpbiB1c2luZzoKCiQkClxiZXRhX3s3fSBccG0gdF97biAtIHB9IFxsZWZ0ICgxIC0gXGFscGhhIC8gMiBccmlnaHQgKSBcY2RvdCBcd2lkZWhhdHtcdGV4dHtWYXJ9KFxoYXR7XGJldGF9X3s3fSl9XnsxLzJ9CiQkCgpgYGB7cn0KaW50ZXIgPC0gcXQoMSAtIDAuMDUgLyAyLCBkZikgKiBzcXJ0KHZhcl83KQpjKGh5cG90aGVzaXNfbW9kZWwkY29lZmZpY2llbnRzWzhdIC0gaW50ZXIsIGh5cG90aGVzaXNfbW9kZWwkY29lZmZpY2llbnRzWzhdICsgaW50ZXIpCmBgYAoKIyMgRXhhbXBsZSAzLjE1IE11bmljaCBSZW50IEluZGV4IOKAlCBQcmVkaWN0aW9uIEludGVydmFscwoKVXNpbmcgdGhlIHByZWRpY3Rpb24gaW50ZXJ2YWw6CgokJAp7eF97MH19XntUfSBcY2RvdCBcaGF0e1xiZXRhfSBccG0gdF97bi1wfSgxIC0gXGFscGhhIC8gMikgXGhhdHtcc2lnbWF9IFxsZWZ0ICggMSArIHt4X3swfX1ee1R9IFxsZWZ0IChYXntUfSBYIFxyaWdodCApXnsxfSB4X3swfSBccmlnaHQgKV57MS8yfQokJAoKV2UgY2FuIG9idGFpbiB0aGUgJFxoYXR7XHNpZ21hfSQgZGlyZWN0bHkgZnJvbSB0aGUgbW9kZWwgb3IgdXNpbmc6CgokJAoobiAtIHApIFxoYXR7XHNpZ21hfV4yID0gXGVwc2lsb25ee1R9IFxlcHNpbG9uCiQkCgpgYGB7cn0Kc2lnbWFfbW9kZWwgPC0gc3VtbWFyeShoeXBvdGhlc2lzX21vZGVsKSRzaWdtYQpzaWdtYV8yIDwtIHQoYXMubWF0cml4KGh5cG90aGVzaXNfbW9kZWwkcmVzaWR1YWxzKSkgJSolIGFzLm1hdHJpeChoeXBvdGhlc2lzX21vZGVsJHJlc2lkdWFscykgLyAobnJvdyhtdW5pY2hfcmVudF9pbmRleF8wMSkgLSA4KQpzaWdtYSA8LSBzcXJ0KHNpZ21hXzIpCmMoc2lnbWFfbW9kZWwsIHNpZ21hKQpgYGAKCldlIGFyZSBhbHNvIGdvaW5nIHRvIHVzZSB0aGUgZGVzaWduIG1hdHJpeDoKCmBgYHtyfQpkZXNpZ25fbWF0cml4IDwtIGFzLm1hdHJpeChoeXBvdGhlc2lzX21vZGVsJG1vZGVsKQpkZXNpZ25fbWF0cml4WywgMV0gPC0gMQpjb2xuYW1lcyhkZXNpZ25fbWF0cml4KVsxXSA8LSAiMSIKaGVhZChkZXNpZ25fbWF0cml4KQpgYGAKCmBgYHtyLCBmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9NiwgZmlnLmFsaWduPSdjZW50ZXInfQphcmVhX3NlcSA8LSBzZXEobWluKG11bmljaF9yZW50X2luZGV4XzAxJGFyZWEpLCBtYXgobXVuaWNoX3JlbnRfaW5kZXhfMDEkYXJlYSksIGxlbmd0aC5vdXQgPSAxMDApCmFyZWFpbnZjIDwtICgxIC8gYXJlYV9zZXEpIC0gbWVhbigxIC8gYXJlYV9zZXEpCgpiZXRhbSA8LSBoeXBvdGhlc2lzX21vZGVsJGNvZWZmaWNpZW50cwp5ZWFyY2MgPC0gbXVuaWNoX3JlbnRfaW5kZXhfMDFbbXVuaWNoX3JlbnRfaW5kZXhfMDEkeWVhcmMgPT0gMTkxOCwgXVsxLCBdCiMgVGhyZWUgbGFzdCB0ZXJtcyB3ZXJlIGFkZGVkIGZvciBjb21wbGV0ZW5lc3MKY29uc3RhbnQgPC0gYmV0YW1bMV0gKyBiZXRhbVszXSAqIHllYXJjYyR5ZWFyY28gKyBiZXRhbVs0XSAqIHllYXJjYyR5ZWFyY28yICsgYmV0YW1bNV0gKiB5ZWFyY2MkeWVhcmNvMyArIGJldGFtWzZdICogMCArIGJldGFtWzddICogMCArIGJldGFtWzhdICogMAoKcGxvdCgKICBtdW5pY2hfcmVudF9pbmRleF8wMSRhcmVhLAogIG11bmljaF9yZW50X2luZGV4XzAxJHJlbnRzcW1fZXVybywKICB4bGFiID0gJ2FyZWEgaW4gc3FtJywKICB5bGFiID0gJ2VzdGltYXRlZCByZW50IHBlciBzcW0nCikKYmV0YV8xIDwtIGh5cG90aGVzaXNfbW9kZWwkY29lZmZpY2llbnRzWzJdCmxpbmVzKGFyZWFfc2VxLCBjb25zdGFudCArIGJldGFfMSAqIGFyZWFpbnZjLCBjb2wgPSAncmVkJywgbHdkID0gMikKCmNvbmZfaW50ZXIgPC0gcXQoMSAtIDAuMDUgLyAyLCBkZikgKiBzcXJ0KGNvdl9tYXRyaXhbMiwgMl0pCmJldGFfY29uZiA8LSBjKGJldGFfMSAtIGNvbmZfaW50ZXIsIGJldGFfMSArIGNvbmZfaW50ZXIpCmxpbmVzKGFyZWFfc2VxLCBjb25zdGFudCArIGJldGFfY29uZlsxXSAqIGFyZWFpbnZjLCBjb2wgPSAnZGFya2JsdWUnLCBsd2QgPSAyKQpsaW5lcyhhcmVhX3NlcSwgY29uc3RhbnQgKyBiZXRhX2NvbmZbMl0gKiBhcmVhaW52YywgY29sID0gJ2RhcmtibHVlJywgbHdkID0gMikKCnhfbyA8LSBhcy5tYXRyaXgoY2JpbmQoCiAgMSwKICBhcmVhaW52YywKICB5ZWFyY2MkeWVhcmNvLAogIHllYXJjYyR5ZWFyY28yLAogIHllYXJjYyR5ZWFyY28zLAogIDAsCiAgMCwKICAwCikpCgpwcmVkX2ludGVyIDwtIE5VTEwKCmZvciAoaSBpbiBzZXFfbGVuKG5yb3coeF9vKSkpCiAgeF9vX2kgPC0gYXMubWF0cml4KHhfb1tpLCBdKQogIHByZWRfaW50ZXIgPC0gYyhwcmVkX2ludGVyLCBxdCgxIC0gMC4wNSAvIDIsIGRmKSAqIHNpZ21hICogc3FydCgKICAgIDEgKyAoCiAgICAgIHQoeF9vX2kpICUqJSBzb2x2ZSh0KGRlc2lnbl9tYXRyaXgpJSolZGVzaWduX21hdHJpeCkgJSolIHhfb19pCiAgICkpCikKCmxpbmVzKGFyZWFfc2VxLCBjb25zdGFudCArIGJldGFfY29uZlsxXSAqIGFyZWFpbnZjIC0gcHJlZF9pbnRlciwgY29sID0gJ2dyYXknLCBsd2QgPSAyKQpsaW5lcyhhcmVhX3NlcSwgY29uc3RhbnQgKyBiZXRhX2NvbmZbMl0gKiBhcmVhaW52YyArIHByZWRfaW50ZXIsIGNvbCA9ICdncmF5JywgbHdkID0gMikKYGBgCgpUaGlzIGNhbiBhbHNvIGJlIGRvbmUgdXNpbmcgYHByZWRpY3RgOgoKYGBge3IsIGZpZy53aWR0aD0xMCwgZmlnLmhlaWdodD02LCBmaWcuYWxpZ249J2NlbnRlcid9CnhfbyA8LSBhcy5kYXRhLmZyYW1lKGNiaW5kKAogIGFyZWFpbnZjLAogIHllYXJjYyR5ZWFyY28sCiAgeWVhcmNjJHllYXJjbzIsCiAgeWVhcmNjJHllYXJjbzMsCiAgMCwKICAwLAogIDAKKSkKY29sbmFtZXMoeF9vKSA8LSBjb2xuYW1lcyhoeXBvdGhlc2lzX21vZGVsJG1vZGVsKVsyOjhdCmNvbmZpZGVuY2VfaW50ZXJ2YWwgPC0gcHJlZGljdChoeXBvdGhlc2lzX21vZGVsLCB4X28sIGludGVydmFsID0gImNvbmZpZGVuY2UiKQoKcGxvdCgKICAgIG11bmljaF9yZW50X2luZGV4XzAxJGFyZWEsCiAgICBtdW5pY2hfcmVudF9pbmRleF8wMSRyZW50c3FtX2V1cm8sCiAgICB4bGFiID0gJ2FyZWEgaW4gc3FtJywKICAgIHlsYWIgPSAnZXN0aW1hdGVkIHJlbnQgcGVyIHNxbScKKQpsaW5lcyhhcmVhX3NlcSwgY29uZmlkZW5jZV9pbnRlcnZhbFssICJmaXQiXSwgY29sID0gJ3JlZCcsIGx3ZCA9IDIpCmxpbmVzKGFyZWFfc2VxLCBjb25maWRlbmNlX2ludGVydmFsWywgImx3ciJdLCBjb2wgPSAnZGFya2JsdWUnLCBsd2QgPSAyKQpsaW5lcyhhcmVhX3NlcSwgY29uZmlkZW5jZV9pbnRlcnZhbFssICJ1cHIiXSwgY29sID0gJ2RhcmtibHVlJywgbHdkID0gMikKCnByZWRpY3Rpb25faW50ZXJ2YWwgPC0gcHJlZGljdChoeXBvdGhlc2lzX21vZGVsLCB4X28sIGludGVydmFsID0gInByZWRpY3Rpb24iKQpsaW5lcyhhcmVhX3NlcSwgcHJlZGljdGlvbl9pbnRlcnZhbFssICJsd3IiXSwgY29sID0gJ2dyYXknLCBsd2QgPSAyKQpsaW5lcyhhcmVhX3NlcSwgcHJlZGljdGlvbl9pbnRlcnZhbFssICJ1cHIiXSwgY29sID0gJ2dyYXknLCBsd2QgPSAyKQpgYGAKCiMjIEV4YW1wbGUgMy4xNiBQb2x5bm9taWFsIFJlZ3Jlc3Npb24KCgo=